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ABSTRACT 


Ongoing interest in establishing a base on Mars has spurred a need for regular and 
repeated visits to the red planet using a cycling shuttle to transport supplies, eq uipment 
and to retrieve surface samples. This thesis presents an approach to designing an optimal 
heliocentric cycling orbit, or cycler, using solar sails. Results show that solar sails can be 


used to significantly reduce V,s at Mars and Earth. For example, using a reasonably 


high performance solar sail, a V_|,,,.= 2.5 km/s is possible at every synodic period using 


Mars 
a two-dimensional orbit model. Lower performance sails were also modeled resulting in 


paths that behaved more like a ballistic Aldrin cycler with higher Vs. Double 


rendezvous missions were explored where the spacecraft must match the velocities of 
both Earth and Mars, offering promising trajectories for Mars sample return missions. 
The solutions to these missions, although not necessarily cyclers, show that using a sail to 
rendezvous with and remain near Mars for an optimal amount of time will minimize the 
total transit time between Earth and Mars. Géeneral-purpose dynamic optimization 
software, DIDO, is used to solve the optimal control problem using a pseudospectral 


method using both two- and three-dimensional elliptic orbit models. 
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I. INTRODUCTION 


Spacecraft orbiting along cycling trajectories between Earth and a nearby planet 
such as Mars could be tremendously valuable for missions ranging from sample returns 
to ferrying supplies and personnel between the planets like escalators. Naturally 
occurring cycler trajectories would be ideal, requiring no added energy other than that 
provided by the gravity fields of the target planets. If such “free ride” paths exist, they 
could host spacecraft making an infinite number of round trips without requiring 
propulsion systems. This hypothetical path would be characterized by a purely natural 
Keplerian trajectory that repeats itself endlessly interrupted only by planetary swingbys. 
Such cyclers could be analyzed for various target planets, however Earth -Mars cyclers 
will be studied in this thesis because of relevance to establishing a base on our 
neighboring red planet. To this end, two prominent Ear th-Mars cycler concepts have 
been proposed- the Aldrin cycler' and VISIT cycler”. The Aldrin cycler uses gravity 
assists from both planets and small well-timed AVs to maintain a continuous cycling 
trajectory. While the revisit times are appealing (7 round trips in 15 years), on-board 
propellant is required to provide a 15 year cumulative AV of approximately 2 km/s. 
Because no fuel is expended trying to reduce speed near a planet, hyperbolic excess 
speeds are high, in excess of 5 km/s. A VISIT cycler, on the other hand, uses neither 
gravity assists nor fuel for AV burns for orbit shaping, but rather resides in a natural 
heliocentric orbit that makes regular passes of both Earth and Mars. The advantage to 
this cycler is that no propellant is required to maintain the orbit for up to 20 years, but the 
primary disadvantages are the rather long revisit times (3 Earth visits and 4 Mars visits in 
15 years) and large approach distances needed to avoid unwanted gravitational 
interaction. The purpose of this paper is to present an Earth-Mars cycler concept using 
solar sails that has the advantages of both cyclers without the on-board propellant. 

Because it has been shown that a nearly ballistic cycler orbit can be maintained 
for 15 years’, it is reasonable to speculate that a solar sail cycler is not only feasible, but 
more capable since it provides free controls. With solar radiation pressure provid ing free 


thrust to the solar sail, a cycler orbit may be maintained with short revisit times as well as 


slow planetary approach speeds. Small approach velocities to Earth and Mars at short 
distances are desirable for operational reasons and because they exhibit attractive orbit - 


shaping characteristics due to large angle swingby deflections. 


Given that a solar sail requires no propellant, we are free to determine a cost 
function that does not minimize the fuel. What then defines an optimal solar sail cycle r? 
Often the time of flight or an undesirable cycler characteristic like large approach 
velocities may be minimized. A design parameter of the sail itself, such as the size of the 
sail may be the set cost function. This thesis presents a framework for solving such 
optimal low-thrust cycler trajectories using the general-purpose dynamic optimization 
software DIDO’. It is worth mentioning that although solar sails are chosen as the 
propulsion method and cycling orbits are selected as the mission for this thesis, the 


framework presented here is applicable to any type of propulsion system and mission. 


Il. APPROACH TO OPTIMAL CONTROL SOLUTIONS 


A. ROADMAP 


At the outset of this thesis work, there were no known works related to optimal 
Earth-Mars solar sail cyclers, therefore a cautious approach to the optimal control 
solution was necessary. Starting with “simple” two point boundary value astrodynamic 
problems, a systematic approach was used that gradually increased the model complexity 
to arrive at solutions. These solutions were used to both validate the numerical 
optimization method and to obtain milestone profiles providing initial guesses and 
bounds to more complex problems. Benchmark optimal control problems were solved 
using DIDO, and the output solutions were compared with known solutions. Benchmark 
problems included Mars flyby, Mars rendezvous and Earth-Mars-Earth double 
rendezvous missions. Numerical results were also propagated using initial conditions and 
the output control histories for further validation. For many of the simpler problems (i.e. 
no intermediate event conditions) a costate profile was obtained corresponding to all the 
Lagrange multipliers at certain points in time along the path. Using the costates with the 
derived solar sail steering control law enabled a check of compatibility between the 
DIDO costate and output control histories. The resulting Hamiltonian output was also 
available to check for optimality. After the benchmark problem solutions had been 
verified, additional constraints and event conditions were included to obtain the more 


complex cycler trajectory that includes planetary swingbys of Earth and Mars. 


The above approach was used for both two-dimensional and three-dimensional 
models. The two-dimensional solutions approximated both Earth and Mars orbits as 
being circular and coplanar serving as an initial estimate for the higher fidelity models. 
Fortunately, the relative inclination of Mars’ orbit with respect to the Earth ecliptic is 
only 1.85° and both orbit eccentricities are not too large, so this approximation is not too 
bad. Increasing the fidelity of the cycler model, Mars was given an elliptic, but coplanar 
orbit to observe the differences in the optimal paths. Finally, Earth and Mars orbits were 
inclined and both eccentricities were taken into account. The results from the battery of 


2-D problems provided bounds and initial guesses for their 3-D counterpart problems. 
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Comparison of the 2-D and 3-D solutions to the same optimal control problem provided 
some interesting insights into astrodynamic optimization. 


1. Coordinate Systems 


To improve the performance of the numerical analysis, it served well to choose a 
coordinate system with states that change slowly with time and dynamics that have few 
or no singularities. Using a coordinate system in which states vary relatively slowly with 
time seemed to allow accurate solutions with fewer optimization “node” points, or points 
along the path that are used to approximate the continuous states and controls in the 
numeric optimization process. Additionally, singularity avoidance is necessary to obtain 


feasible solutions. 





Figure 1 Elliptic Orbit in Polar Coordinates 


For the two-dimensional model, a heliocentric polar coordinate system was used 
as shown in Figure 1. Because the maximum solar sail thrust is inversely proportional to 
the distance from the sun, this coordinate system was well suited for this problem. The 
only singularity in the polar equations of motion occurs at the center of the sun - an easily 
avoidable situation with proper state bounds (see Appendix A). Velocity components 
were expressed in terms of radial velocity and along-track velocity (ref 4 p. 43). The 
along-track velocity will be referred to as the transverse velocity and is defined as being 
normal to radial position vector in the orbital plane. Since the target orbits were 
approximated as circular and coplanar in the 2-D model, there was no need to define a 


principle direction (i.e. along vernal equinox). 
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Likewise, for the 3D model spherical coordinates were chosen to represent 3-D 
states and dynamics because the resulting optimal state and control profiles were 
expected to vary slowly with time allowing excellent numerical approximations using 
relatively few nodes. Furthermore, these coordinates made it very simple to model the 
thrust due to the sail since thrust is dependent on radial distance, r (see Appendix B). By 
comparison, states of orbiting bodies in Cartesian coordinates oscillate. This could have 
lead to an undesirable condition in case there turned out to be too many cycles in the 
solution for the number of nodes to approximate accurately. It should be kept in mind, 
however, that Cartesian coordinates present perhaps the best choice when more 
perturbations are desired in the dynamics model (in this case, more nodes should be 
used). Equinoctial elements provide a set of singularity-free slowly changing variables, 
but because it was desirable to keep boundary conditions simple and “intuitive”, they 


were not used (ref 4 p. 142). 


2. Scaling 


With the coordinate system chosen to reduce the number of optimization nodes 
required to accurately approximate states and controls, it became imperative to select 
units that would reduce numerical conditioning by avoiding the extremely large and small 
numbers. There is no need to confine the states to familiar units such as miles, 
kilometers, seconds, etc. A useful unit system for modeling interplanetary travel is the 
canonical system. Choosing a distance unit (DU) to be 1 A.U. and the normalized solar 
gravitational parameter u to be 1, other state units may be calibrated accordingly. All 
states and controls have units that are of the same order of magnitude and dimensio ns are 
scaled by a factor of some characteristic dimension. This non-dimensionalization process 
also enhances visualizing metric relationships. For instance, the circular Earth orbit 
radius is 1 DU while the Mars orbit radius is 1.524 DUs, over 50% greater than Earth’s. 


Heliocentric canonical units are summarized in Table 1. 


Dimension Non-dimensional unit | mks equivalent Comments 


Distance unit DU 1.496E8 km Semi-major axis of Earth orbit = 1 
AU. 
circular orbit at 1 AU 


VU = DU/TU 29.78 km/s Velocity of body in 1 AU circular 
orbit. 
Gravitational u = 1 DU/TU 1.327E11 km’/s Solar gravitation parameter 
parameter 





Table 1 Canonical Units Used to Non-Dimensionalize the States and 
Parameters. 


3. Bounding the Problems 


Before approaching the full cycler model accounting for gravity assists, simpler 
two-point boundary value problems were analyzed. As already mentioned, this permits 
validation of the optimization software. Additionally, important insights may be gained 
by observing the behavior of simpler benchmark problems. For both the 2-D and 3-D 
models, three benchmark problems were solved prior to the coding of the cycler models. 
These are the Mars flyby, Mars rendezvous, and Earth -Mars-Earth double rendezvous 
round trip missions. The Mars flyby mission is a simple minimum time trajectory from 
Earth orbit to Mars orbit. This solution provides the absolute shortest time that a solar 
sail of a given performance can reach the Mars orbit. Other missions that include an 
Earth to Mars transit may use this time as a lower bound since it is impossible to reach 
Mars in less time. A Mars rendezvous mission is one in which the solar sail is tasked to 
reach Mars at its local orbital velocity (V..= 0 with no Mars gravity acting) in minimum 
time. The resulting profile is not only useful in observing optimal rendezvous behavior, 
but also for bounding the double rendezvous problem that followed. The double 
rendezvous mission required the spacecraft to rendezvous with Mars, then track the 
planet for a time (with no effect due to Mars gravity), and then return to rendezvous with 
Earth. This mission introduced phasing constraints to ensure that Earth was there when 
the spacecraft returned. With the double rendezvous bounds established in this manner, a 


state discontinuity was introduced representing a gravity assist from the Mars encounter 
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making a zero-sphere patched conic cycler model complete. At this point, all the models 
and conditions had been established for a solar sail cycler and coded for use with the 


numeric optimization software, DIDO. 


4. Numerical Analysis, terminology and DIDO* 


Ideally, the solution to the optimal trajectory would view the entire space-time 
path of a spacecraft at once and simultaneously locate the optimal position for each point 
along that path meeting certain boundary conditions and dynamic constraints (Figure 2 
and Figure 3). Limiting ourselves to current computer methods, however, we settle for 
using certain points along the space time path that will enable us to approximate the 
continuous path, thus dscretizing the trajectory into a workable form for a non-linear 
programming (NLP) solver. An Euler method would provide an approximation that 
would be burdened by large errors. It was desired for this complex problem to use the 


best numerical approximations possible in solving the optimal control problem (OCP). 
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Figure 2 Path and Target manifolds (orbit positions) 


space-time plot 


time [yrs] 





Figure 3 Path Meeting Event Manifolds (Planet Encounters) in Space (x-y) and 
Time. 


The numerical solution to the optimal control problem presented in this thesis 
uses the software DIDO to interface with a generic third party NLP solver. The NLP 
solver requires an NLP variable vector, constraints and a cost function. DIDO provides 
these requirements set up in a numerically optimal fashion. All state and control path 
information are captured at unevenly distributed node points, thus the entire dynamic 
problem can be viewed as a static problem with the complete trajectory history at selected 
node points manifested as a single variable vector. The node points are selected in such a 
manner that the state and control approximation error is minimized. These numerically 
optimal node points are called Legendre-Gauss-Labatto (LGL) points and correspond to 
the zeros of the first derivative of the Nth degree Legendre polynomial>*®. The 
continuous states and controls are approximated by using the variables at the 
predetermined LGL node points and Lagrange interpolating polynomials, which are well 


suited for interpolation between unevenly distributed data points. 
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Dynamic constraints are imposed by forcing the derivatives of the approximated 
states to equal user-defined differential equations (equations of motion) at all the node 
points. Obtaining the derivative of the state variable vector is accomplished using a 
differentiation matrix defined in ref 5. This process of discretization and differentiation 
is called the Legendre pseudospectral method and has been used successfully in a variety 


of applications (refs 5-8). 


Initial, intermediate and final event conditions ae supplied by the user and 
structured into a boundary condition vector that DIDO provides as constraints to a third - 
party NLP solver. The approach used in this numerical solution to optimal control does 
not require propagation or shooting type methods. An initial guess is required, however 
it does not need to be a good one, or even a feasible one. Typically the end conditions 


are supplied and DIDO interpolates between these points as a guess. 


To assist in the reading of the numerical methods used in this thesis, some 


definitions from ref 3 are supplied below. 


Definitions: 


Nodes are discrete points along the path where states and controls are defined. 
The entire continuous path is represented by the states and controls at these discrete 
points. The states and controls at the all the node points form part of a vector making up 
the NLP variable that inputs to the NLP solver. Node locations are unevenly distributed 
and are pre-determined at the LGL points. LGL points occur at the zeros of the first 
derivative of the Legendre polynomials in order to minimize the mean -squared error in 


Lagrange interpolating polynomial approximation of states and controls. 


Knots are double node points that are part of the optimization process. Trajectory 
end points are also called knots. Interior knots are used to define intermediate events 
where there may or may not be state discontinuities. LGL node distribution always 
concentrates more nodes next to knots. Intermediate or interior knots have a left and 


right side where state constraints are imposed. 


Events are the set of conditions that define the manifolds which constrain the path 


(at the knots). The event manifolds usually consist of the boundary conditions; in this 
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thesis they represent the locations and velocities of Earth and Mars along their respective 
orbits in space and time. Within the DIDO structure, knots are the primary way to 

represent an event condition. An example of event boundary conditions is shown below 
where e represents the event vector, and e; and e, represent the lower and upper bounds of 


the vector respectively. 


[te | [to | Oa | 
Vo Vor Vou 
e=|r, | where e, =| r, | and e, =| r;, 
Vr Ve Vex 




















such that 


e, Se(x(f,),x(t),X(t,),p t,t, t, Se, 


where the vector p contains any static parameters that may constrain the events. Events 


form the constraints that DIDO supplies to the generic NLP solver and are listed in tables 


for the various OCPs in this thesis. 


Path constraints, g , are the state-control constraints that couple control and state 


variables in the following manner. 
8, <g(x(1),Wt),1)<g, 
For numerical reasons it is sometimes advantageous to represent the solar sail 
controls as the sine and cosine components of the cone and clock control angles instead 
of the actual angles themselves. The path function in this case may be represented as 


2 


2 2 

uy tu, —1 : . ‘ 

g -| peas | where g, =g, =0 forming an equality constraint. The actual control 
ux, +uy — 


angles may be recovered by taking the arctangent of coupled controls (using Matlab’s 


atan2 function serves well here). 
Bounds are the upper and lower limits on states, controls or times. The resulting 
restriction is in the form of an equality or inequality constraint. For instance, a set of 


state variables may be bounded by xe S CR” usually written in the form xek , x,,] 


or x, <x<x,,. Equality constraints are implemented simply by making x, =x 
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Within DIDO, states and controls are discretized, so variables are bounded at the node 


points. 


Event Constraints are those linear or non-linear conditions that restrict states at 


the events. For instance, in an orbital rendezvous problem, the final spacecraft velocity is 


coupled with the final radius of the target planet according to the relation vy, ices 
Vr 
forming a non-linear event constraint. When a state at a particular event condition is 
unconstrained, it is considered free. For example, when intercepting Mars in minimum 
time, we do not specify where in space or time the intercept must occur, thus final time 


and final angular displacement are free variables. 


B. VALIDATING SOLUTIONS 


Once a solution to a particular trajectory was obtained, it was important to see if 
the path was physically feasible and if so was it optimal. Checking for feasibility proved 
to be relatively straightforward, while verifying optimality was a bit thornier. This 
process linked the stark numerical optimization solution with a physical spacecraft 


trajectory that “made sense”. 


1. Feasibility 


A glance at the state data at the knots would verify that the event constraints had 
been met. In order to verify that the entire state history conformed to physical laws 
governed by the equations of motion we use a propagator. Taking the control history 
generated by the numeric optimization and the initial conditions, the path of the sail could 
be propagated by means of a numeric ordinary differential equation solver on the non- 
linear equations of motion. The propagation used the same dynamics employed by the 
DIDO code to verify that DIDO was properly applying the dynamic constraints to the 
path. To verify the equations of motion in the dynamic model used by DIDO, the 
equations were propagated without controls to generate familiar Keplarian orbits. The 


state and control history figures in this thesis generally show the DIDO states represented 
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by markers at the discrete node points and the propagated solution as a line that will pass 


through DIDO’s markers when the solutions are consistent. 


2. Optimality 


After testing for feasibility, the solution was checked for optimality. The primary 
method of verifying that an optimal trajectory had been achieved was observing the 
behavior of the Hamiltonian, H(x, u, 4, t). DIDO-derived Hamiltonian values were 
generated at the predetermined node points for simple problems (i.e. no intermediate 
knots) and were used to verify minimum time optimal controls. Since none of the 
problems posed in this thesis have a Hamiltonian as explicit function of time, then we can 


say 


Hamiltonians that are constant with respect to time are necessary (but not 
sufficient) conditions for optimality in these cases. Moreover, the optimal control 
solution is the one where H is minimized with respect to controls, u thus defining an NLP 
problem with the Lagrangian of the Hamiltonian (for details see Appendix A). For 
minimum time problems, H =—1 for all time. The first and second order necessary (but 
not sufficient) conditions for optimality using the Hamiltonian are written as 

OH _ o°H 


—=0 and 


ou ow es 





Additionally, dual variable outputs were used with the derived solar sail control 
law to produce controls that could be compared with the DIDO -derived controls. This 
would test if the optimal controls were conforming to the analytical optimal control 
steering law solution. Occasionally there existed one or more local minimums in which 


case varying the initial guess assisted in finding a better local minimum. 


3. Comparison with Other Optimal Solutions 


Often it serves well to verify DIDO-generated optimal solutions that use the 
pseudos pectral method with solutions created using other methods. As mentioned before, 
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there are no known optimal control solutions to the actual solar sail cycler, however there 
are solutions to the simpler “benchmark” problems. One of the reasons for using the 
methodical build-up approach using simple building blocks in the form of benchmark 
problems was to verify that the coded models were producing verifiable results before 
proceeding to the uncharted waters. Comparisons are presented in the benchmark 


problems when other solutions are available. 
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Il. DEVELOPING THE OPTIMAL TWO-DIMENSIONAL 
CYCLER 


A. TWO-DIMENSIONAL MODELS 


1. Sail Model 

A solar sail is a very large thin lightweight structure that transfers momentum 
from inbound solar photons to the spacecraft. Newton’s laws dictate that the sail will 
accelerate twice as fast when photons are reflected than when they are absorbed. Thus, 
for maximum effectiveness the sail must be made of a very reflective material. For this 
analysis, the sail is modeled as a flat film (non-billowing sail) with a perfectly specular - 
reflecting surface on one side only. To characterize sail performance, we use the sail 


lightness number, B , defined as the ratio of the acceleration from solar radiation pressure 


(SRP) to the acceleration from the sun’s gravity (ref 9 pp. 38-40). 


2W.R. 
: Zee Acos’o 
B = srp = Mc 
A oray i 


where W, is the solar energy flux at 1 AU, R, is the Sun-Earth distance, A is the sail 
area, O is the angle of the sail with respect to the Sun and w is the solar gravitational 
parameter’. Both accelerations are proportional to the inverse square of the radial 
distance from the sun, so 8 is an apt indicator of sail performance independent of 
location. A lightness number of $B =.17 was used to model a reasonably high 
performance solar sail for most problems in this thesis for comparison to other solutions 
using the same sail!°. Often the characteristic acceleration!', d,, is used to describe the 


performance of the sail instead the sail lightness number. This characteristic acceleration 
is the change in velocity that the sail would experience at a one AU distance from the sun 
when the sail exposes all of its area toward the sun. The relationship between the 


lightness number and the characteristic acceleration is 
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pe =a, 


% 


For example, a sail lightness number off =.17 corresponds to a characteristic 


. mm 
acceleration of a, = 1—. 


S 


The sail control angle, ©, is defined as the angle between the sun-sail position 
vector, r, and the sail normal, n called the cone angle (Figure 4). This angle determines 
both the magnitude and direction of a solar force imparted to the spacecraft. Notice that 
when & =0, maximum thrust is achieved and when & = 90’, thrust is zero. The sail can 
never impart a force toward the sun (only gravity can do this). To ensure the sail will not 


flip over exposing its non-reflective side, the control angle is bounded by 


T T 
1 -—s<sas— 
(1) - 5 
r 
n 
O 
Figure 4 Sail normal, solar radial vector and control angle o. 


2. Dynamic Model 
Using the polar coordinates in the heliocentric frame, the two-dimensional 
equations of motion are derived in Appendix B Given the state vector x=[r 9 uv], 


the equations of motion are expressed as x =f (x uw) where control u=Q, 
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The states r, 8, u, and v are the radial distance, angular displacement, radial 
velocity and transverse velocity respectively. The angular displacement, 0, is decoupled 
from the other equations and may be excluded in special cases making xe R*. For 
numerical analysis of some problems, © was retained as a state for proper intercept 


phasing during encounter events (more on that later). 


Having thus imposed the dynamic constraints between the desired manifolds, we 


seek to constrain the conditions at the manifolds. This is accomplished by setting the 


constraints of the desired events. 


3. Events Model 


All the benchmark missions modeled here share the same cost function, dynamic 
equations of motion and sail control model. The distinguishing feature of each mission is 


its particular event function. The events are arranged into an event function vector e such 
that ee R”’ where N, is the number of event constraints. Elements of e may be linear 


or nonlinear expressions and may be bounded as equality or inequality constraints. These 
event constraints for the flyby and rendezvous missions are summarized in the following 


table where a, is the semi-major axis of Mars’ orbit (1.524 AU). 
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Table 2 Event Conditions for Flyby and Rendezvous Problems 





Notice that all the initial states (denoted by the ‘0’ subscript) are equality 
constraints with the same lower and upper bounds ( //b,ub]) indicating that the spacecraft 
will launch from a circular Earth orbit at Earth’s orbital speed. The final rendezvous 
along-track velocity (denoted with the “f’ subscript) is the only non-linear constraint. 
These events are assembled within DIDO into a _ constraint — vector, 


e€=[7, Dy Uy M% Ty 9 pu ive Te —1]’. The last element is written in a different form to 


avoid the +/- ambiguity using the square root. 


B. THE OPTIMAL CONTROL PROBLEM 


We now come to posing the solar sail optimal control problem (OCP). This is 
where one must question what exactly are we trying to do. With most conventional 
propulsion schemes, it is desirable to minimize fuel consumption. Solar sails, however, 
require no propellant so the goal of most of the solar sail trajectories presented her e is to 


accomplish a given mission in minimum time. The general OCP is stated as follows’. 


Minimize the Bolza cost functional: 
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Subject to: 
Dynamic constraints x =f (x(t),u@) 
where f is given in equation (2) and U are the controls, 
Path constraints g, < g(x(t),u(2),4) Sg, 


Event constraints e€, Se(X,, X,,X, 7, ¢,¢ p p) <e, 


and bounds on state and control variables 


X, S x(t) <x, 


u, Su(t) <u, 


With an established coordinate system and scaling system, we are postured to 
model the sail, dynamics and events. These models are structured into a format ready for 
DIDO to use as the OCP constraints. Before attempting a numerical solution, we desire 


the analytical solar sail steering law to help verify the solution as described below. 


C. SOLAR SAIL CONTROL LAW 


Following the guidelines set forth by the minimum principle (Appendix A), we 
start by establishing the cost function. To minimize time, we choose to express the cost 


functional from the previous section in Lagrange form. 


E=0,F =1 


The state space vector and its derivative are written as 


and f(x,u)= Bs 
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where the control is given by the sail cone angle a(r). 


Notice that 8 is uncoupled from the other equations, therefore it will be ignored 


for the rest of the sail control analysis. Thus, we write the Hamiltonian as 


2 
H(x,u,A) =F +A'f =1+Au+A, Ep Ben'a ra [Ht B ost sine 
r r ror 


where A is the costate vector and the solar gravitational parameter has been normalized 


to unity. 


The Lagrangian of the Hamiltonian is therefore 


2 
L(x,u,A)=H+ug=1l+Aut+, = Sr Bcosta |, [- B cos’ sino } [Ol 
r r r 


7 7 

where Hg, is now the covector associated with the path constraint “55 as . A 
Applying the minimum principle yields 

= =-A, Ecos’ sina +A, [Boose sin? + F cos‘ I LL, =0 


SO Hy may be written as 


B 


Ll, = cosa Ba cosa sina +20 sin’ a —A, cos’ a ) 
r- u v v 


For the given control constraints, the KKT conditions give 


Limiting the analysis to the interior controls where = <a< . and U,=0, we can 
obtain the control law. 
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(3A,, cosa sin & + 2A, sin? —2, cos*a) = 0 


Dividing the equation by cos’a produces a quadratic equation with tana. as the 


independent variable. 


2A, tan? + 3A, tana -A, =0 


a‘ + Jon, +82, 


» An, 


Vy 





This tangent control law will not have quadrant ambiguity since we limited the 
control angle by =5S as .. 


D. RESULTS 


Before analyzing results, it is useful to correlate certain key sail angles with 
physical meaning. As mentioned in the discussion of the sail model, the thrust magnitude 
due to the solar radiation pressure (SRP) on the sail is not a free control variable as it is 
with many conventional propulsion systems. The acceleration imparted by the sail, T ,/m, 
has a magnitude that is dependent on radial distance and the control angle, a, according 
to the relationship 2 = BE cose . Notice that for a given r, maximum radial thrust 
occurs when = 0 and it’s sail area is totally exposed to the sun. Alternatively, the sail is 
effectively “off’ and the trajectory becomes ballistic when a = +7/2. The sail can never 
isolate tangential thrust from radial thrust. Radial thrust is present for all control angles 
except for & = +7/2 when the sail edge is toward the sun. To assist in interpreting the 
control profile, we determine the control angle at which maximum transverse acceleration 


occurs. The transverse acceleration is given by: 


Bu 


Les. IS Tt Tt 

a, =— sina =—cos"a sina, —-—<a<— 
2 

m r 2 2 


Setting the derivative with respect to @ to zero for a given r yields the control 


angle providing maximum transverse acceleration. 
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da : 
—t=~—2cosa sin*a+cos*a =0 
(4) do. 


O =.615= 35.2° 


1. Benchmark Problem Solutions with Circular Coplanar Orbits 


a. Earth-Mars Flyby 


The purpose of solving the flyby mission serves several purposes. First, 
the solution may be compared to other known solutions to verify that the numerical 
analysis using DIDO works properly. Second, the minimum time of flight serves as a 
lower time bound for future coplanar trajectory optimization problems. Third, 
groundwork is established to model a bounded initial C3 from Earth allowing a limited 


“boost” from a conventional rocket at time t = 0. 


Given an initial guess of the initial and final positions and velocities, 
DIDO outputs the minimum time Mars flyby states and controls (plotted as markers in 
Figure 5). For a minimum time flyby path between circular coplanar orbits, a sail with 
lightness number B=.17 takes 0.45 years. The probe sails past Mars at a speedy 8.7 km/s 
(relative to Mars). Sail attitude favors a large local transverse acceleration profile 
initially to rapidly build radial acceleration, and then gradually exposes more sail area to 
the sun throughout the maneuver until the sail normal is parallel with the sun’s rays 


(Figure 6). DIDO costate outputs representing 4 =[A,, ,.A,,. a] are plotted in Figure 7, 


and may be used with the tangent steering control law in equation (3) to generate derived 
controls. The true anomaly costate, A, , is zero as expected since @ does not appear in 


the Hamiltonian. A comparison of costate-derived controls with the DIDO output 


controls appears in Figure 8. 
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Mars min time flyby. DIDO vs. propagated states and controls. 16 nodes. 
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Figure 5 States and Control for Mars Flyby Mission. Markers are DIDO 
output and lines represent propagated path. 
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Figure 6 Mars Flyby Trajectory with Sail Profile and Normal. 


Costates for Mars minimum time flyby problem. 16 nodes. 
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Figure 7 Mars Flyby Mission Costates 
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Figure 8 Mars Flyby Cone Angle Control 


Thus far, the flyby model has restricted the sail to start under its own 
power without the aid of a kick motor. Modeling a kick motor is accomplished by 
changing the initial conditions in Table 2 for the u and v velocity components to be 


constrained variables according to the following relationship. 


O<IV,/< AV, 


where V, =[u,.v,]'-—V,. Not surprisingly, the optimal solution makes use of whatever 


AV... 18 permitted to intercept Mars in the quickest time, although the direction of 


max 


departure depends on the maximum allowable kick. The mission is accomplished much 


faster than the previous sail-only solution. For example, when AY, ag the optimal 
s 


solution uses all the initial help it can get reducing the time to intercept to only .22 years. 
We find, however, that the same behavior does not apply to the rendezvou s problem with 
its final velocity constraint. 

Feasibility of the solution is demonstrated in Figure 9 depicting the path 
propagated using an ODE solver given the initial conditions and the DIDO -derived 
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control profile. Although difficult to rigorously prove that this is the optimal trajectory, 
we can at least show that necessary conditions for optimality are satisfied in accordance 


with the discussion in the section on Optimality. The Hamiltonian is constant with 


2 
respect to time (Figure 10) and the second order condition 





= 0 is true (Figure 11). 


7s 





Figure 9 Mars Flyby Mission Propagation (line) with DIDO Output (dots) 
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Hamiltonian Profile for Mars Flyby Mission 
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Second Order Condition for Mars Flyby 
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b. Rendezvous 


As expected, the Mars rendezvous mission takes longer to accomplish 
than the flyby mission, requiring a full 1.11 years, which is consistent with the solution 
presented in ref 12, (Figure 12). The spacecraft final velocity is constrained to match 
Mars’ orbital velocity, which costs the sail extra transit time to accomplish. A surprising 
feature of the trajectory is that the probe sweeps an arc outside of the circular orbit of 
Mars (Figure 13). This maneuver is observed in more complex minimum time solar sail 
missions as well since it appears to make the best use of radial and transverse thrust. 
There is a point on the outbound path (almost a quarter of the way into the mission time) 
when the control angle rotates from a small angle, where much of the sail area is exposed 
to the sun, to a more edge-on aspect. The sail attitude gradually rotates after this to favor 
more transverse acceleration as it sweeps past Mars’ orbit radius and approaches Mars. 
By the time the sail has reached Mars, the control angle is 5° greater than its maximum 


transverse velocity setting. 


EM Rendezvous with 50 nodes 


normalized units 
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Figure 12 DIDO State and Control Output (Markers) and Propagated Path 
(line) for EM rendezvous. 
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EM rendezvous with 50 nodes 


ae 





Figure 13 Sail Trajectory and Profile for Earth-Mars Rendezvous 


Interpreting the costate history provides a more complete physical analysis 
of the optimal trajectory. Each costate in Figure 14 represents a Lagrange multiplier that 
signifies instantaneous sensitivities of the cost function to instantaneous variations in the 
corresponding state. Once again, the costate corresponding to 9 is zero since this state is 
completely uncoupled from the other states in the dynamic equations and never appears 
in the Hamiltonian. Several observations may be made regarding the critical point in the 
first quarter of the trajectory when the control angle makes a radical change. Recalling 


the optimal sail steering law in equation (3) it is seen that when i, =0, then the 


V2 


maximum transverse velocity profile occurs. This value would occur when tana = a 


(a =35.2°). Also, when A, >0 then & 90°. Observing the costates for the Earth- 
Mars (EM) rendezvous, we note that A, crosses zero at t = 1.225 years and soon 
afterwards i, reaches its closest point to zero. This is consistent with a rapid change of 


sail attitude at that time during the mission and corresponds to when the sail reaches 


maximum radial velocity. 
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Lagrange multipliers along with state outputs are fed into the derived 
tangent steering law (equation (3)) to obtain a control profile. The resulting control 


profile is shown in Figure 15 and is plotted with the DIDO control output for comparison. 


Costates for EM rendezvous mission. 50 nodes. 
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Figure 14 Costates for Earth-Mars Rendezvous Mission 
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Control Profile 


ido output 
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Figure 15 Comparison of DIDO Controls and Tangent Steering Control Law for 
EM Rendezvous. 


Propagating the spacecraft with the given control history produces the path 


shown in Figure 17 with DIDO states at node points shown for comparison. Testing the 





2 
second order necessary conditions for optimality we observe in Figure 16 that : = 20 
] 


for all time and also the Hamiltonian is fairly constant (Figure 18). 
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Figure 16 Second Order Necessary Condition : z 
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Figure 17 Propagated Path (line) with DIDO State Output (dots) for EM 
Rendezvous 
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50 nodes. 


(ND units) 
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Hamiltonian 
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Figure 18 Hamiltonian for Minimum Time EM rendezvous. 


What happens when we reverse the problem such that we seek the 
minimum time from Mars (V.lmars = 0) to Earth rendezvous? It turns out that swappin g 
initial conditions with final conditions generates a time-reversed state history and a 
control profile that is inverted and reversed. The experiment may be carried out within 


the DIDO structure by simply reversing the event conditions. 
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Table 3 Mars-Earth Rendezvous Event Conditions. 


The state history of the optimal Mars-Earth rendezvous is a mirror image 
of the Earth-Mars rendezvous (Figure 19 and Figure 20). Time of flight, At, and the 
change in angular displacement, A®, are precisely the same for both profiles. 
Furthermore, the control profile is inverted and time reversed. This result leads to an 
even more interesting question; what happens when we desire the sail to launch from 
Earth, rendezvous with Mars, then immediately start its trek back home to rendezvous 


with Earth? This benchmark mission is the subject of the next section. 
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Solar sail trajectory with sail normal vector history 





Figure 19 Mars to Earth Rendezvous. Shows reverse trajectory of EM 
rendezvous mission 
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Figure 20 ME Rendezvous States and Controls 
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C. Double Rendezvous 


As a next step toward modeling a cycler orbit, we explore the double 
rendezvous. In this mission we seek the minimum time that it takes to make an Earth - 
Mars-Earth (EME) round trip with a zero relative velocity at both planets for a given sail 
lightness number 8. To study the behavior of the optimal trajectory, one assumes that 
Earth or Mars gravity will not significantly influence the craft. Low relative velocity is 
desirable while encountering the target planets since it requires less energy for a greeting 


taxi craft to intercept and dock with the passing spacecraft. 


For the preliminary analysis, the initial, intermediate and final orbits are 
modeled as circular and coplanar. Earth phasing for the final Earth rendezvous event is 
accomplished by coupling the sail and Earth final positions using the Earth’s mean 


motion and the time of flight. 


The solution to the minimum time EM rendezvous problem provides a 
lower bound for the time that it must take to reach Mars during the double rendezvous 
problem. Since the fastest possible time to rendezvous with Mars with the given solar 
sail (B=.17) was determined to be 1.11 years from the previous section, this serves as the 
minimum time to encounter Mars halfway through the double rendezvous mission. 
Remembering that the reverse Mars-Earth (ME) rendezvous took exactly the same 
amount of time, we now have bounded the time to complete the whole EME double 
rendezvous. Setting the intermediate destination at the Mars orbit and the final 
destination with Earth and constraining the event velocities to match those of the 
respective planets we are in a position to solve the optimal EME double rendezvous. 
Shown in Table 4 is a summary of the event constraints where superscripts indicate 
before (-) or after (+) the intermediate knot representing the Mars encounter. There is a 
single non-linear event constraint that is responsible to ensure proper phasing (i.e. Earth 
makes more orbits around the sun than the sail does, but still meets with it in the end). 
The output states and controls are shown Figure 21 with the path shown in Figure 22. 
Essentially it turns out to be the EM minimum time rendezvous solution patched together 
with the ME rendezvous solution with slightly longer transit times. The most obvious 


features of the profile are the state symmetries and control antisymmetry about the mid - 
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maneuver point. Event conditions at the interior knot required only that the states be 


continuous and equal to Mars’ state. 


| Event | Mbub] | ub] [lb,ub] [lb,ub] | Event | | Mbyub] ub] 
Pepa ss re fica [or 


Non-linear cons raint cos, +1, (¢; —t)—3,) [1,1] 





Table 4 Event Conditions for EME Double Rendezvous 


For the given sail, the total time to complete the trajectory is 2.41 years 
with both legs of the journey taking 1.205 years. The reason that the EME problem takes 
longer than the patched EM-ME problem is that the final Earth rendezvous event must be 
phased with Earth’s position at the final time. The individual EM and ME rendezvous 
solutions only served to provide a lower time bound for the EME problem, so imposing 
the constraint that the sail final position is phased with Earth is about 10% more than the 


lower bound. 
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EME rendezvous - min time 


—t— radial dist 
—* radial vel 
—+ transverse vel 
—t. sail angle 





Time in years 


Figure 21 States and Controls for EME Double Rendezvous in Minimum Time 


(B=.17) 





Figure 22 EME Double Rendezvous Trajectory and Sail Profile (B=.17) 
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This round trip solution only takes into account a quick trip out to graze 
the Mars orbit and then a return trip home. It may be desirable to have the probe linger 
with Mars for a span of time allowing for some sort of mission while traveling with the 
planet then returning back to Earth. The next question one may ask is the following. 
Given a sail of certain performance and Mars stay time, what is the minim um transit time 
EME trajectory (not including the Mars stay time)? For a given Mars stay time, some of 
the transit time to and from Mars is expended in performing a phasing maneuver to meet 
Earth at the final time. This characteristic phasing maneuver is exhibited in all double 
rendezvous problems with varying stay times but one — the one where the stay time spent 
with Mars happens to be just the right time to allow a return trip to Earth with no phasing 
maneuver required at all. As mentioned, this minimum possible return trip time is 1.1 
years for sail lightness number B=.17. The following analysis shows that given the 
optimal time and angular displacement required for the EM rendezvous (and thus the ME 
rendezvous), we can deduce the Mars stay time that minimizes transit time. The optimal 
transit time and true anomaly traversed by the sail for an EM rendezvous (and ME 


rendezvous) are given by the following. 
t =1.116 years 6° = 4.337 radians 


The following equations capture the motion of the Earth relative to that of 


the cycler craft where t,= Mars stay time, AO,= change in Earth angular position and 


n,and n, represent Earth and Mars mean motion respectively. 
(5) AO, =n,(2f +t,) 


Equation (5) indicates that Earth changes position by A®, during the same 


time that the sail transits outbound in f° years, stays with Mars for ¢,, years, then returns 


home in f years. Recognizing the n, = Pe we can rewrite equation (5) as 


(6) AO, -t,, = 20° 
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During the round trip mission time the spacecraft traverses 0° radians on 


the outbound leg, n,,t,, radians along with Mars, and @° on the inbound leg. 


(7) AO =20°+n,t,, 


All the while, Earth traverses the same span plus an extra N revolutions. 
(8) AG, =AO+20N, (N, =[0,1,2...]) 


Constraining the solution to include only the first Earth pass on the return (N=1), we 


substitute equation (7) into equation (8) and rearrange to obtain 
(9) AO, —n,,t,,= 20°+ 20 
Recognizing that the scaled mean motions of Earth and Mars are given by 


1 
n, =1 and n,, =,/—+ =.534 respectively, we may solve for A@, and ¢,, simultaneously 
a 


using equations (6) and (9). This results in a Mars stay time of ¢, =7.0114TUs = 117 
days and an Earth span of A®,=16.03 radians corresponding to 2.55 years of total 


mission time. 
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Figure 23 EME Double Rendezvous with Stay Time at Mars 


When using the outbound and inbound transit times as a cost function, 
DIDO produces the output shown in Figure 23. The path mimics a patched minimum 
time EM-ME solution where the optimal Mars stay time is 123 days for the given sail 
(5% difference from the estimated value of 117 days). Any more or less time spent at the 
planet will either require a time-wasting phasing maneuver to allow Earth to catch up for 
a rendezvous, or a maneuver to catch up to the speeding target Earth. A plot showing the 
impact of various Mars stay times on mission transit times appears in Figure 24. Notice 
that it is safer to design a shorter-than-optimal Mars stay time into the mission. In the 
event of a schedule slide, the return transit time looks increasingly grim beyond the 


optimal point because the sail has a lot of space to cover as it attempts to catch up to 


Earth. 
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round trip transit duration for various Mars visit times 
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Figure 24 EME Double Rendezvous Transit times for Various Mars Stay Times. 
(B=.17) 
The analysis thus far has only dealt with sails of lightness number B=.17. 
An interesting question that may be asked now is the following. Given a desir ed time of 
round trip flight with no stay time at Mars, what is the minimum size solar sail required? 


This question will be addressed later. 


2. Benchmark Problem Solutions with an Elliptic Coplanar Orbit 


Running the same battery of problems from the previous section using a higher 
fidelity model provides enormous insights into the characteristics of an optimal 
trajectory. Modeling trajectories between circular Earth and coplanar elliptic Mars orbits 
generates results that may be compared with the circular orbits model to reveal the “knees 
in the curve” that the optimization process seeks in an effort to reduce the cost function 
just a little more. Knowing these characteristics of optimal solar sail trajectories will 


assist in understanding the behavior of optimal solar sail cyclers. 


Because the final target manifold represents Mars’ elliptical orbit, the new events 


model must relate the final Mars radial position, r, os and velocity, V_, with the final sail 


m? 


state, X,. 
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Constraining the initial and final sail angular positions to be coincident with the 
initial Earth and final Mars angular positions respectively, we obtain the following 
relationship for coplanar orbits. 

6, =9 ,-(@, +Q,) 
where 8, is the final sail position, @,, is Mars’ argument of periapsis and Q,, is the right 


ascension of the ascending nodes. 


In a perifocal system, Mars final polar coordinates r,,and 8,,, are related by 


Pin 
10 =—P» __ 
mo) cae et €,, C089, 


where p,,in the numerator is the “parameter” or semi-latus rectum of the Mars orbit. 
The final event condition for radial position of the sail is simply 


(11) 1, =Ny 


For the rendezvous mission, we need to target the final Mars velocity as well. The speed 


of Mars at the final time is expressed as 


where a,, is the semi-major axis of Mars. Components of the planetary velocity vector 


are defined in terms of the flight path angle B,, by 





and sinB,, = et sin (0, ) 


cos B,,, = 
Ting Y mf ming 
where h,, represents the magnitude of the Mars angular momentum vector. 


The velocity vector resolved in the local vertical local horizontal frame (LVLH) 


produces 


Unp = Ving sin B.,, 
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Vir = Vip COS be. r 
. . oye _ T _ 
Representing final Mars and sail velocities by V, fp = ay Ving | and V a= [u Vy 
respectively, we obtain an additional constraint for the elliptic rendezvous mission 


(12) V,=V,, 


It now remains to run the same battery of battery of benchmark problems with the 


elliptic Mars event model. 


a. Earth-Mars Flyby 


Modeling Mars’ orbit as an ellipse using eccentricity e, =.0935 provides 


target radial distances that vary with the planets true anomaly (equation (10)). Earth’s 
orbit eccentricity is only .0167, so the circular Earth orbit assumption is a good one 
(elliptic Earth orbit is considered in the 3-D model in the next chapter). Not surprisingly 
the time optimal mission selects a path that drives the sail from Earth orbit to Mars 
periapsis at high speed (Figure 25). This optimal path exploits Mars’ eccentric orbit and 
presents the shortest distance from Earth to Mars, and thus the shortest transit time 
required since the final velocity is unconstrained. The time to intercept Mars using the 


standard sail (8 =.17) is only 137 days, about 27 days faster than the corresponding 


circular orbits model (Figure 26). 


16 nodes. 





Figure 25 Minimum Time EM flyby with Mars Elliptical Orbit. Final position 
is at Mars periapsis. 


EM minimum time flyby with elliptic coplanar orbits. 16 nodes. 
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Figure 26 EM Flyby with Mars Elliptic Orbit 
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b. Rendezvous 


To achieve the rendezvous mission in minimum time the path must 
account for target velocity as well as distance, both of which vary with respect to target 
planet true anomaly. Initial and final true anomalies, thus launch windows and arrival 
dates, are optimization parameters that seek to reduce the cost function by finding the 
optimal points on the target manifolds to bound the path. It is fascinating to observe that 
the optimal rendezvous with Mars (elliptic orbit) is actually faster than the simplified 
circular coplanar rendezvous (0.977 years vs. 1.11 years). Comparison of this state 
history with the circular coplanar counterpart gives insights as to why this occurs ( Figure 
27). In the circular coplanar model, any mission start date is as good as any other 


therefore the boundary condition 8, =0 was a valid restriction where the Mars lead angle 


was determined in the numerical solution to the OCP. Because in the new events model 


we allowed the optimal 6, to be determined in the OCP solution, a launch window was 


chosen such that Mars would be near the slowest point in its elliptic path at the time of 
rendezvous. Mars matches the sail velocity at the final event as it slows down on its 
approach to aphelion. The optimal rendezvous trajectory takes advantage of Mars’ slow 
velocity as the sail approaches the planet at a radial distance that is 95% of the aphelion 


distance. Figure 28 shows how the sail meets the Mars position on its orbit. 


Another interesting departure of the elliptic orbit solution from the circular 
solution is that the spacecraft does not follow a trajectory that sweeps out to a maximum 
radius and then returns inward to meet Mars at the required velocity (Figure 27). Mars’ 
orbital path however does cross inside the spacecraft path as Figure 28 readily reveals. 
The optimal 3-D trajectory takes the fastest path from the optimal location on Earth’s 
orbit to intercept Mars without having to undergo a negative radial velocity during the 


whole maneuver, i.e. u > 0 for all time (Figure 29). 


46 


Optimal state paths for elliptic and circular orbits 
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Figure 27 Comparison of State Profiles Using the Circular (markered lines) and 


Elliptic (thick lines) Orbits Models. 
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Figure 28 EM Rendezvous with Mars Elliptic Orbit 
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normalized units 
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Figure 29 EM Rendezvous with Coplanar Elliptic Orbits. 


Cc. Double Rendezvous 


Once again, aphelion provides the optimal Mars rendezvous point (Figure 
30 and Figure 1). The total transit time is only 2.3 years in contrast to the 2.4 years it 
took to reach the circular Mars orbit. In contrast to the double rendezvous orbits (see 
Figure 21), the radial velocity, “, does not change as much in the vicinity of the interior 


knot thus “flattening” the radial distance profile as the spacecraft sweeps past Mars. 
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EME double rendezvous between elliptic orbits. 50 nodes. 
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Figure 30 States and Control for EME Double Rendevous with Elliptic Mars 
Orbit. 
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Figure 31 EME Double Rendezvous Trajectory and Sail Profile with Elliptic 
Mars Orbit. 
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3. Earth-Mars Synodic Cycler Solution 


Having established a collection of solutions for basic solar sail trajectories, we are 
in a position to take the next step toward a cycler model. There are two key differences 
between the double rendezvous mission of the previous section and the solar sail cycler 
mission. With cyclers, end conditions are equal (i.e. cyclic) by definition to ensure 
repeatability. Additionally, gravity assists at target planets will be modeled. Planetary 
swingbys offer enhanced cycler performance since the turn angles are optimized to create 
a path that achieves the minimum cost. Since the time to complete a cycle is 
predetermined by the synodic period (igure 32), a minimum time synodic cycler 
problem has no meaning for spacecraft shuttling between circular coplanar orbits. This is 
not so for the cycler between non-circular coplanar orbits. A new cost function is desired 
that is not burdened with minimizing fuel or time. Because an attractive feature of 
cycling trajectories is a slow swingby velocity at Earth and Mars, these velocities formed 


the cost function. 


a. Events Model 


The key to modeling a repeating cyclic trajectory is to constrain the initial 
state of the spacecraft to equal the final state of the spacecraft. Since the Earth and Mars 
orbits are approximated as circular and coplanar, the problem is simplified in that the 


initial relative angular position between Earth and Mars (lead angle ¢ ) may be used to 


constrain final angular position for a single cycle since planetary distances and velocities 
are independent of the their inertial angular positions (Figure 32). Lead angle is 


determined from the event s tates as follows. 
O. =0, 5 +7, (2 —85) 
8. =6,+n,¢,-%) 
where 


8.0 =0; 4 n,,(t; = ty) 
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Mars and Earth angular (8, and @,) are obtained using their respective 


m 


mean motions (n,, and n, ). 
The Earth-Mars lead angle at initial and final times are 


Ce = 8, —9 
CG =8,7 -9. 


Initial and final relative angular positions are constrained by 
(13) cos(6, -¢,) =I 
where C, =€,-2nN, , N, =0,1,2.... 
Numerically, equation (13) is preferable over the above equation since the 
cosine function will permit multiple revolutions without introducing integer variables. 
The cyclic end condition for the spacecraft radial distance is expressed 
simply as 
(14) ty =r 
Initial and final velocities are also constrained to be cyclic and will be 
addressed later. 
Given the condition in equation (13), the final time, t,, depends only on 


the relative mean motions of the planets given by the synodic period, 


2n 








For Earth and Mars, the synodic period is 2.135 years'*. Although not 
explicitly constrained to the synodic period in the numerical analysis, the resulting final 
time must equal this synodic period in order for the cycler to be periodic. Having 
established cyclic end conditions, we now turn to setting up gravity assist conditions for 


planetary encounter events. 
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Figure 32. Two-Dimensional Earth-Mars Cycler Geometry with Circular 
Coplanar Planetary Orbits 

Force imparted by the solar sail shapes the path between event manifolds, 
but gravity assist maneuvers at the event manifolds drive the form of the whole cycler 
trajectory. As with the conventional Aldrin cycler, gravity assists are implemented in the 
solar sail cycler to shape the trajectory at these planetary encounter events to improve the 
revisit times. To model the swingby events, state discontinuities are employed at the 
planet to change the velocity direction and magnitude of the spacecraft in the heliocentric 


frame (sometimes called the “zero sphere patched conic’’* 


or matched asymptotes 
model). It is assumed that the interaction time with the subject planet is negligibly small 
in comparison to the total cycle time. The velocity of the spacecraft with respect to the 


planet changes direction such that V,,* = V,, +AV, where V,, =[u,v]'—V, and V, is 


the velocity of the subject planet relative to the sun. The position states, r and 0 are 
constrained to be continuous at both Earth (equation (14)) and Mars encounter events 


given below. 


(5) 
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The AVs due to the swingbys are optimally chosen in the OCP solution, 
but they must be properly constrained. Constraints on the velocity changes are most 
easily imposed using the planet frames where inbound and outbound velocity magnitudes 


are equal just before and after a planet encounter forming the event condition, 


(16) lV, i] = ve =V,. (planet frame) 





where V, is the hyperbolic excess speed. The velocity direction change is expressed in 


the planetary frame using the turn angle, 5 , which exists in the region shown in . 
Velocities before and after a swingby event in the planet frame are coupled by the cosine 


of the turn angle, 


- + 2 
V,, *V,, =V., cosd 


- + 
e 


V, Sp 
(17) cos6 = ee (planet frame) 


co 


The spacecraft experiences a direction change during the interaction that is 
restricted by the hyperbolic excess speed and the permissible periapsis pass distance from 


the center of the subject planet, r,. This restriction is expressed by the following 


relationship in the planet frame where V,, is scaled by the circular orbit speed at the 


surface of the planet and r, is scaled by the radius of the planet (ref 15 p. 24). 


: 1 
(18) sin (5 ye eS ae ee 
Although equations (17) and (18) themselves are not event conditions, 


they limit the achievable change in velocity direction, AV, due to the swingbys at Mars 


and Earth. Design limitations include a selected r,,,,, that is well above the atmosphere 
where drag effects are negligible and an r,... such that the encounter occurs close enough 
to the planet to execute a desired task. Fora given V_, the turn angle is maximum when 


ro=r 


Pp pmin 


and minimum when r, =r,, 


ma: 


.- Substituting these values into equation (18) and 


solving for 6 as a function of V_, provides an expression for 6, and 6. ,, respectively. 


a 
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To include the case of maximum deflection in the opposite direction it is necessary to 


place the lower bound of acceptable turn angles at —-d_,.and -6,,. (Figure 33). With 


max 
equation (17) characterizing the instantaneous change in path direction and equation (18) 
providing the limits, the boundary conditions are expressed at Earth and Mars event 


manifolds as 


(19) cos(6.,,,,) <cosd <cos(5,,,,) 


cos(8) feasible region 





Figure 33 Gravity Assist Geometry. 


Velocity constraints at Mars and Earth have identical form, however the 
initial velocity magnitude at Earth has an additional limitation based on available 
departure rocket capability. Because the sail’s journey starts at Earth, the initial 
conditions are bounded by maximum C; available. Presumably an impulse rocket is used 
to start the solar sail craft on its cycler trajectory, so the initial velocity relative to Earth, 


V, | 


0 ‘earth? 


is restricted. The magnitude of V,|,,,,, 18 limited by a maximum allowable 


earth 


velocity change at t=O provided by a kick motor, which provides another boundary 


condition 


(20) 0<|V, | 





<AV (Earth frame) 


max 
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earth 


The direction of V, | 


learn 18 driven by the optimal control problem and is 


limited only by the allowable turn angle at each Earth swingby. 

Finally, phasing the spacecraft with Earth at the end of a cycle was 
considered. Earth encounter events were constrained to ensure that the sail trajectory 
intersected Earth’s orbital path at precisely the time that the planet is at that same 
location. The circular orbit assumption is particularly useful for ensuring proper phasing 
of events since the angular position of a planet is a linear function of time. The final 


angular position is given by 
(21) 6, =8,, +1, -20N,, 


where N, is an integer number of Earth orbits, n, is Earth’s mean motion, and @,, is the 
angular position of Earth at t=0 (when we the Earth position is coincident with the sail, 
0,, =9,). Having established the sail, dynamic and events models, we set up the optimal 


control problem. 


b. Solar Sail Cycler Problem Formulation 


The solar sail cycler optimal control problem is constructed using 
weighted spacecraft V.s at the Mars and Earth encounters as the cost function. The 
optimal path is subject to two-dimensional equations of motion, cyclic end conditions, 
and planetary gravity assist constraints. The cost function uses a parameter y to weight 
the V..s while the initial and final state and control conditions are constrained to be equal 
to ensure repeatability of the trajectory. We can write the optimal control problem as the 


following. 


Minimize the cost 


(22) J(x;, Y) = YV.. Lae +(1 =H) Vi. ees 


Subject to dynamic constraints from the equations of motion (equation (2)) 


and event constraints that model repeatability (equations (13) and (14)), continuity 
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(equation (15)), swingby effects (equations (16) and (19)), launch limitations (equation 


(20)) and phasing (equation (21)). 


In addition to the bounds on controls (equation (1)), there are bounds on 


states. States are bounded only to avoid singularities in the dynamics. 


For this paper, a single Earth-Mars cycle was modeled where the 
following initial, intermediate and final conditions in astronomical units define the target 
manifolds 

mH =r, =1 
r, =1.524 
The parameters that bound the path deflections at the event manifolds are 


the maximum G available at the initial Earth orbit departure, AV. 


max 


in canonical units, 


and the minimum and maximum periapsis pass distances at Earth and Mars. 


AV ax, = 2 (~6 km/s) 
lamin = 1-06 Mars radii 


r = unbounded 


mmax 


r,. =1.16 Earth radii 


emin 


r,_. =10 Earth radii 


emax 


Cc. Synodic Cycler Results and Analysis 


Shown in Figure 34 is the state and control angle output from DIDO for a 
single cycle of the synodic cycler with y=1. As expected, the time required to complete a 
cycle under this set of constraints was 2.135 years, the Earth-Mars synodic period. A 
quick glance at the state history reveals that the spacecraft sails from 1 AU out to Mars’ 
orbit at 1.524 AU and then returns to Earth. A discontinuity in both velocity states occurs 
at 1.524 AU and 1 AU representing Mars and Earth swingbys respectively. The optimal 
path makes use of large gravity assist maneuvers (in the planet frame) during planet 
encounters owing to slow V..s. The spacecraft initially accelerates in the radial direction 
while decelerating in the transverse direction. At the appropriate time, the sail rotates to 
an attitude that favors more and more positive transverse acceleration to intercept Mars 


with the lowest V..to minimize the cost function while setting up for the swingby event 
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initiating the return trip. Following Mars swingby, the spacecraft sweeps out to nearly 2 
AU to ensure proper phasing for Earth intercept. Sail attitude gradually reaches a 
maximum negative transverse acceleration profile @ = -35°; see equation (4)), then 
“shuts off’ and follows a ballistic path as it presents an edge aspect to the sun. A plot of 
the sail trajectory with Earth 1, Mars 2 and Earth 3 encounters is shown in Figure 35. A 
similar gravity assist is accomplished at the Earth encounter and, since cyclic end 
conditions were imposed, the same control profile will reproduce the trajectory 
repeatedly. Owing to these constrained end conditions, the initial Earth departure 
hyperbolic excess velocity only required 4.3 km/s — not the maximum allowable limit of 


6 km/s. 


Synodic Cycler with 115 total nodes 


ND units 
cone angle [deg] 





Time [years] 


Figure 34 DIDO States (markers) and Control with Propagated Path (line 
through markers) for a Single Synodic Cycle. 
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Solar sail trajectory with sail profile 
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Figure 35 Single Cycle Path of Solar Sail Cycler with Minimum V..|nars 


A noticeable difference between the solar sail cycler and the traditional 
impulse rocket Aldrin cycler is the large swingby angular deflections with respect to 
Mars and Earth. The Aldrin cycler, because it is minimizing fuel mass, resides in a 
natural Keplarian orbit most of the time. As such, it tends to have a large V.. in excess of 
6 km/s at Mars and excess of 5 km/s at Earth, thus restricting turn angles. The solar sail, 
on the other hand, can change orbital energy with no impact to the cost function and 
achieve low hyperbolic excess speeds that permit large turn angles. The results of this 
analysis show that a 75° Mars turn angle and a 29° Earth turn angle provide the optimum 
path. Interestingly, the sail never goes down to the minimum allowable pass distance 
with Mars to get a bigger swingby deflection. Furthermore, the sail swings by Earth at 
the maximum allowable perigee distance, not the minimum allowable distance. Table 5 


summarizes the cycler parameter data for the cost function with y=1. 
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Table 5 Earth-Mars Solar Sail Cycler Flyby Data (6=.17, y=1) 





It should be noted that the parameters shown in Table 5 are highly 
sensitive to the states right before the Mars and Earth encounters. These numbers 
represent the optimal parameters given the approximated state and control history. By 
slightly modifying the approximations using more node points, states preceding a knot 


could change enough to generate a slightly different optimal parameter set. 


To verify the solutions, states were propagated using initial conditions, 
gravity assist conditions and the DIDO-generated control history using a Matlab® ODE 
solver as described in the Validating Solutions section. Propagated states are shown 


passing through the DIDO output markers in Figure 34 producing the path in Figure 36. 


ND 





Figure 36 Propagated Path (line) with DIDO State Output (dots) for Single 
Cycle. Solution uses 115 total nodes (45 before and 70 after the interior knot). 
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To capture the effectiveness of the cyclic end condition constraint, the 
initial and final propagated states were compared for a single cycle. A comparison 
summary appears in Table 6 where the error between final and initial conditions of the 


propagated path are shown with states given as radial distance, r, Mars lead angle, C , 


velocity magnitude, V , and zenith angle, y (complement of flight path angle). 


92.38° 93.86° 1.60 


Table 6 Cyclic End Condition Errorsfor Synodic Cycler 





Up to this point, we have only minimized the hyperbolic excess speed at 


Mars with the weighting factor y in equation (22) equal to unity. To suit the needs of 
any particular cycler mission, however, it may be desirable to minimize a combination of 


at both Earth and at Mars. Varying the weighting factor produces the range of V..s 


shown in the graph in Figure 37. A y of approximately 0.3 will minimize the cost 


function the most with least total V... 


V-infinity at Earth and Mars 
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Figure 37 Varying Y in the Complex Combination Cost Function 
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d. Remarks 
Using a reasonably high performance solar sail to achieve optimal 
heliocentric cycler trajectories is a viable alternative to traditional impulse rocket cyclers. 


There are several interesting ways to pose an “optimal” solar sail cycler problem. One 


such problem would be a synodic cycler that achieves Mars and Earth encounters with 


the constraint that V.=0. This double rendezvous synodic cycler could also pose an 
intriguing design optimization problem in which one desires the minimum sail lightness 
number required to achieve a double rendezvous in a synodic period. These cases are 


investigated in the next section. 


4. Fun with Cycler Trajectories 


a. Double Rendezvous Synodic Cycler 
Using the standard performance sail (B =.17), a cycler has been 


presented that minimizes V.. during planetary encounters. It may be interesting to seek a 
synodic cycler that is constrained to rendezvous with each planet such that V.. is zero at 
both planets. Recall that performing an EME double rendezvous mission with the 
standard sail with unrestricted end conditions took at least 2.41 years. This sail would 
never reach Earth again within a synodic period of 2.135 years. A higher performance 
sail however might be able to. Intuition tells us that there ought to be a sail with just 
barely enough area to achieve an EME double rendezvous within a synodic period. This 
sail design optimization problem may be stated as follows. 


Find the minimum sail lightness number that would enable a double 


rendezvous cycler where V.. = 0 at both Earth and Mars. 
Minimize: J=8 
Subject to dynamic constraints 


x-f(x,u)=0 
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and event constraints given in Table 4 with the additional constraint that the end 
conditions are cyclic. Hyperbolic excess velocity is zero at the initial and final times, so 


the only additional constraint is that the EM lead angle, C , is the same at both times. 
me = C f 


The solar sail lightness number is established as a static parameter, a sail 
design characteristic that remains constant in time. The results in Figure 38 and Figure 
39 show that the path has characteristics of both a synodic cycler and a double 
rendezvous. There is symmetry in the r and v states, and antisymmetry in the u state and 
control angle like the double rendezvous. Initial and final conditions are precisely the 
same as in a cycler, although no gravity assists are used. The time of flight turns out to 
be, of course, the EM synodic period. The minimum sail lightness number required to 
perform such a mission is B=.297. This corresponds to a sail that is nearly double the 
size of our standard solar sail with the same payload mass! 


Minimum Sail Lightness Number for EME Double Rendezvous 


ND units 





Time [years] 


Figure 38 Minimum f Solar Sail States and Controls for an EME Double 
Rendezvous 
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Solar sail trajectory with sail profile history. 
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Lightness no = .3 (min beta to achieve EME double rendezvous) 


Figure 39 —‘ Trajectory for Minimum B EME Double Rendezvous 


b. Taxi Propellant Cost 


One mission that could benefit from a cycling orbit is replenishing 
supplies at a station on or around Mars. In concept, a taxi craft could leave its parking 
orbit about Mars and greet the passing cycler sail on its hyperbolic trajectory around 
Mars. A better cost function in this case would be the fuel required to meet the cycler on 
its swingby path from a parking orbit (equivalently, we could minimize the Av). We 
assume a circular parking orbit with a radius equal to the closest cycler approach 
distance. Since the same exact pass conditions would be met every cycle, it is presumed 
that the taxi craft is positioned in this orbit. Following the discussion in ref 16 p. 102, the 


Av required to go from a circular orbit to match the cycler’s hyperbolic orbit is 


=, pm 
pm 


where r,,, andv,,, are the periapse position and velocity with respect to Mars. Since the 


energy of a hyperbolic orbit is 
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Enyp =o = 


hone 


and the periapse velocity of the spacecraft with respect to Mars is 


The change in velocity required to join the cycler craft on its hyperbolic 


| 2u,  [ 
= 2 
Vom [Vee oF See fe 
r r 
pm pm 


Note that if the circular parking orbit grazes the closest approach point of 


trajectory is then 


the cycler orbit (for fixed r,) the Av is proportional to v,. Thus for this taxi model, a 


pm 
minimum v,,, is equivalent to a minimum taxi Av problem. For elliptic or 3D target 
orbit models where the closest cycler pass distances vary from visit to visit, the optimal 
taxi intercept path would differ. In this case, the cost would also have to include 
propellant expended to get from a nominal parking orbit to a non-grazing hyperbolic 


cycler path. 


c. Profiles Using Different Sail Performances 


As with any other blossoming technology, advances in sail material design 
are related to the amount of interest and thus funding. In order to make a cycler mission 
feasible it is useful to know what mission designs are available for any given sail 
performance. In the next analysis, a range of sail lightness numbers were fed to the 
dynamics model in equation (2) and analyzed in the same fashion pres ented using the 


standard sail. Results for y =1 in the cost function of equation (22) for a range of 
lightness numbers are shown in Figure 40. The time to Mars, ¢,, and V_|,_,, are used to 


compare trajectory characteristics. The higher performance sails tend to take a longer 


time to reach Mars in order to reduce the hyperbolic exc ess speed at Mars. 
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Optimal Cyclers with Different Sail Lightness Numbers 
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Figure 40 The Effect of Varying Sail Performances on a Cycler with y = 1 in the 
Cost Function 
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As the sail performance is reduced down to £ =.01, the cycler flyby 
characteristics look more like a ballistic Aldrin cycler that passes closer to Earth using 
more bend angle. The sail attitude only changes the orbit enough to prevent the sail from 


dipping below the minimum pass distance restriction. A sail with 1/10 " the performance 


of the standard sail (i.e., B=.017) will reach Mars in a short time, get a small bend, then 
sweep outside Mars orbital path to return to Earth close enough to get a large bend. The 


planet encounter parameters are shown in Table 7. 


i V..(km/s) | rp (planet radii) 8 | 





Table 7 Solar Sail Cycler Characteristics for Cycler (B=.017, y=1). 
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Time to Mars 





E. VALIDATION OF SOLUTIONS 


The primary method of validation was accomplished through numerical 
propagation. Using the same dynamic equations and initial conditions employed by 
DIDO along with the control history, a propagator will generate a trajectory that may be 
compared with the output state history. A match of resulting states indicates that the 
solar sail with the given control history follows a feasible path conforming to physical 
laws. The control profile may be obtained from the direct DIDO output, or derived from 
the costate history (for OCPs without an interior knot) and the sail control law. The latter 


generally provides “smoother” control for interior controls (when control angle is not at a 
limit, ie. 8 a<%). A comparison of DIDO and costate-derived controls were 


shown in Figure 8for the flyby and Figure 15 for the rendezvous problem. 


Numerical propagation was accomplished using the Matlab ® ODE45 and ODE23 
solvers. Controls had to be interpolated between DIDO node points in order to produce 
an accurate sail attitude at time steps generated by the ODE solver. Generally, a cubic 
interpolation served well while a spline method proved inaccurate in regions with 


concentrated node points, i.e. near knots. 


Validations of the benchmark problems are shown in Figure 9 and Figure 17. 
DIDO paths and the propagated paths closely match as shown in Table 8 listing the mean 
squared error for the different runs (DIDO dynamic constraints were met with less than 


10° for both runs). 


Fit 





1.0952e-009 3.4343e-006 
1.0258e-009 3.4093e-006 


Table 8 Mean Square Error of DIDO States Compared with Propagated 
States Using Matlab® ODE45 
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Propagation of swingby trajectories, such as those used in the cycler problems, 
became more of an art than a science. The propagation was enacted for one cycle only 
with an intermediate event condition defined at Mars’ orbital radius. To simulate the 
gravity assist, a switch was used in the propagator to add a AV to the sail’s velocity prior 
to the encounter equal to the AV determined by DIDO. Because of the sharp 
discontinuity in the path velocity due to the gravity assist at the knot, the subsequent 
trajectory was driven by where the assist occurred. Small differences in where the 
outbound path encountered Mars caused the DIDO states and propagated results to differ 
on the inbound leg Figure 41. The discrepancy lay in how the sail was being controlled 
just prior to the interior knot. Output controls often appear “shaky” near a knot. Wit hin 
the propagator, these somewhat erratic controls must be approximated at the time steps. 
If the time steps of the propagator were large in comparison to the DIDO controls at the 
LGL points, approximation errors resulted. These small errors are enough to slightly 
change the Mars swingby event location causing the remainder of the trajectory to 
deviate from the DIDO solution. The discrepancy is resolved by forcing the propagator 
step sizes to be approximately the same size as the distance between DIDO -generated 
LGL points near the knot (Figure 42). In this way, better approximations are made near 
the defined points producing closer matches between DIDO and propagated paths. This 
can be accomplished by either imposing a maximum step size on the propagator or by 
adjusting the number of nodes used by DIDO to manage node spacing. Better control 
approximations in the neighborhood of the planetary encounters yielded a closer match 
between the optimal solution states and the propagated states. This, however, only 
confirms that the output really does match a propagated solution for the given erratic 
control behavior just before the hard knot. We desire a solution that does not require the 


sail to perform radical attitude changes near the Mars encounter. 
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Propagator 
Mars Encounter 






DIDO Mars —~~~ | 
Encounter 


Figure 41 Propagated Path with “Missed” Mars Swingby Phasing. Propagated 
states use controls interpolated at time steps different than DIDOs LGL distributed 
time steps causing differences at the interior hard knot. 





Figure 42 DIDO Control Output at LGL Node Points and Interpolated Controls 
used in the ODE45 Propagation Near a Knot. Step sizes match fairly well. 


68 


The controls may be made “smoother” by adding “inertia” to the controls; i.e. by 


limiting the rates. The control then becomes the time rate of change of the cone angle. 
where the control is u, =Q(f). 


Using this control and state variable convention provided inherently stable cone 
angle control, a more desirable design for a solar sail attitude control system. 
Additionally, the pitch rate, &, was be restricted to provide the “inertia” to the sail 
without having to switch from a 3 DOF problem to a 4 DOF problem. Numerical 
solution time was expected to shorten as well since the hodograph of “control” is 


convex (see Appendix E). 
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IV. DEVELOPING THE OPTIMAL THREE-DIMENSIONAL 
CYCLER 


We now extend the boundary conditions and dynamics of the 2-D solar sail 
trajectory models of the previous chapter to include a third dimension and use the same 
methodical approach to solve the optimal cycler. The definition of a “cycle” for the 3 
dimensional problem however becomes more complex than the 2-D case. Initial and 
final cycle end conditions do not repeat exactly in a synodic period as they did in the 2-D 
circular coplanar model. Planetary orbit inclination and eccentricity make conditio ns 
such that Earth and Mars do not repeat their relative positions with each other but about 


every 15 years. 


The motivation for obtaining a 3-D optimal control profile for a solar sail cycler is 
not only to increase the fidelity of the model but also to compare results to the 2-D case 
to learn more about the nature of optimal cycler orbits. The inclination of Mars’ orbit 
with respect to the ecliptic plane is roughly 1.85°, not a great difference from the 
coplanar case. Elliptic orbits for Earth and Mars are modeled, so optimal trajectories are 
similar to the elliptic coplanar solutions. We perform the same series of benchmark 3D 


solutions to rate them against their 2D circular and elliptic counterparts. 


A. THREE-DIMENSIONAL MODELS 


1. Sail Model 


Our sail in the new model now has an extra degree of freedom due to the addition 
of the third dimension. With most conventional engines, the controls generally have 
three degrees of freedom, one for each of three dimensions. However, with the sail, 
thrust magnitude is dependent on the cone angle providing a constraint, therefore only 
two degrees of freedom are required. The cone angle serves well as one of the control 
variables since the thrust magnitude is related to it. Another angle is required to 
determine where on the cone the sail normal vector lies. It is convenient to define a clock 


angle as the angle between the projection of the sail normal onto a plane normal to the 
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sun-line and a reference direction in the same plane; so we adopt a model used in ref 9 
p.115. This reference direction is taken to be the vector normal to the instantaneous 

orbital plane of the sail. Figure 43 makes the representation more clear with unit vector 
Pp as the reference direction. Note that a positive clock angle rotation is in the negative r 


direction using the right hand rule. 






Plane perpendicular tor 


a = cone angle 
d= clock angle 


Figure 43 Solar Sail Control Model for 3 Dimensional Dynamics 


2. Dynamics Model 


Since the thrust magnitude is dependent on the radial distance from the sun, 
spherical equations of motion provide a simple way of including the acceleration due to 


the sail (see Coordinate Systems Section) 


For the three degree of freedom state variable x=[r9  v,,@,,@,]’ the state 


derivative is given by 
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—| -2v@ +H cos? cisinacdosd —@,° sing cosd 
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where the last three terms are from the equations of motion (Appendix C). 


In choosing these dynamic equations, we must bound the states to avoid 
singularities at r=0 and cosd =0. It turns out that avoiding these state values keeps the 


Jacobian of f free of singularities as well. 


3. Events Model 


The 3D target orbit manifolds are the locations and corresponding velocities of 
Earth and Mars along their respective inclined elliptic orbits. These events define the 
boundary conditions of the optimal control problem. In defining the boundary 
conditions, it is necessary to know the relative orbital shapes and orientations of the 


departure and destination planets. The orbital elements of Earth and Mars orbits are 


summarized in the following table '*. 


(a ee 





isea | 085 ose | 2865 


Table 9 Earth and Mars Orbital Parameters 
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In all the problems presented, the planet’s true anomaly at the events remains a 


free variable to be determined in the numerical solution. 


In order to rendezvous with a planet, it is essential to describe the position and 
velocity of both spacecraft and planet in a common 3-D (Np=3) coordinate system. 


Although the spacecraft states are represented in spherical coordinates, the events at 1, 
and t, are given in perifocal coordinates in the respective planetary orbital planes. To 


match the sail states with the planetary states at the end conditions, Earth and Mars 
orbital planes were transformed into a common heliocentric -ecliptic coordinate system 
along with their respective velocities along the orbital paths. For simplicity, the frame is 
referred to as the E-frame. These orbital states define the target manifolds to which the 


spacecraft event conditions are to be constrained (see Appendix D). 


First, the initial Earth angular position and final Mars angular position are 
constrained to be equal to the sail’s position in the common E-frame. Resolving these 
manifolds in the Eframe we obtain the planetary positions in a common coordinate 


system. 


Constraining the initial and final sail angular positions to be coincident with the 
initial Earth and final Mars angular positions respectively, we obtain for small 


inclinations (ref 4 p. 135) in the E-frame 


8. =8, -@, +Q,) 6, 9 ,-@,, +2) 
where @p is the initial spacecraft position and @, is the final sail position. 


In a perifocal system, the planetary polar coordinates r, and @, are related by 


Pe Pin 


=-——— r. = 
1+e,cos8,, 7 


To 
1+e, cos0 
m mf 


where p in the numerator is the “parameter” or semi-latus rectum of the respective orbit. 


Now the initial and final sail target positions in Cartesian coordinates are given by 
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ro COS@y) rg COS(O ;) 
R,, =A,| 7,,sin@,) | for Earth and R.,,,, =A,,| 7,, sin(®,) | for Mars 
0 0 


where A, and A, are the respective Earth and Mars 3-1-3 rotations transforming into 


the E-frame (Appendix D). In like manner the velocities of the planets may be defined in 


their own perifocal frame by 





where V,, is the speed of the planet, 





r,| is the distance of the planet from the sun, and a, 


is the semi-major axis of the planet. Components of the planetary velocity vector are 


defined in terms of the flight path angle B by 








cosB,, = Fe and sin, ey: 


| P pop 
where h, represents the magnitude of the planet’s orbital angular momentum vector. 


A planet’s velocity vector resolved in the local vertical local horizontal frame 


(LVLH) has components 
u,=V, sin B., 
v, =V,cos B, 


Finally resolving these components in Cartesian coordinates and transforming 
them into the E-frame yields the target velocities V.z and Vy (transformation is given in 
Appendix D). For the rendezvous mission, Vy," is required to match the final sail 


velocity, but knowledge of Mars’ velocity is not needed for the flyby mission. 


The constraints may now be set such that the spacecraft states at the initial and 
final events are equal to the target planetary states. The initial and final sail positions are 


given as follows. 
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(2 %r | 
R,. =| ¥ Vee 


The left column represents the sail position at the beginning, r,, and the right 
column represents the sail position at the final time, r,. All coordinates are resolved in 


the Eframe. After applying the transformation given in equation (45), the spacecraft 


velocities at events in LVLH frame are given by 


The velocities need to be transformed from LVLH into Cartesian coordinates in 


the E-frame using a transformation matrix B defined as 


B=[B] =[R,(-9)][R,)] 
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where the R matrices correspond to standard rotation matrices (Appendix D). The 


spacecraft velocity in the Cartesian E-frame is therefore 
V.=BV, V,eR% 
The left column of V_, represents the initial sail velocity, V,, and the right 
column of V,,, represents the final sail velocity, V,. Having defined all the necessary 


states in the same coordinate sys tem, it is now a simple matter to set boundary conditions 


for the initial and final rendezvous events. 


The sail position starts at Earth and ends at Mars 


aw 
as) 


and likewise the initial sail velocity is the orbital velocity of Earth and final velocity is 


equal to the orbital velocity of Mars. 


(26) V,- V,,=0 


(27) V,-V,,=0 


The final velocity event constraint in equation (27) applies to the rendezvous 
mission only where the spacecraft must match the velocity of Mars in its orbit. Table 10 
summarizes the event conditions in equations (24) through (27) for a rendezvous. All 
coordinates are resolved in Cartesian coordinates in the E-frame. Within the DIDO 
framework we code the event conditions in a more compact form that ensures end 


condition manifolds are constrained. 
R,,-R,, =0 and V—- Ve, =9 


For only two end conditions without intermediate conditions the dimensions of all 


the matrices are N,, xN,,, where N,, is the number of dimensions and JN, is the number 


of events. 


Flyby constraints | Rendezvous constraints 
Event 
[lower,upper] [lower,upper] 


R,] [Ry Rel 


[R 


LV.->Vi-1 [V5 Ved 


Table 10 Event Conditions for 3-D Flyby and Rendezvous Missions. 
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B. THE OPTIMAL CONTROL PROBLEM 


In keeping with the step-by-step approach to reaching a cycler model, we first 
seek the minimum time solutions to the 3-D flyby and rendezvous missions. To this end 
the new 3D OCP is stated in a similar fashion as the 2D OCP. Because the state variable 


has now expanded to six variables (3 position, 3 velocity), we constrain the state 


77 


derivatives to the 3D equations of motion (equation (23)). States are bounded to avoid 


the additional singularity at the solar poles. 


C. SOLAR SAIL CONTROL LAW 


Since the sail attitude control system is responsible for driving both the cone 
angle, o&, and clock angle, 5, we seek two steering laws. As with the 2-D model, we turn 


to the principles of optimal control theory to generate these steering laws. Tailoring the 


Bolza cost function in Lagrange form for a minimum time problem we get 


rae (Fé 


Li) 


E=0,F =1 


Using states consistent with a spherical coordinate system, 
x=[7,0 ,0 ,v,,@, ,@ ; |’, we construct our state derivative vector f from the equations of 


motion (Appendix C). 


Bu 


7,” COS” > + TW, te a oO 
(28) £(x 02,8) =|» 8 9 2 


| 2 (ro, sing —v, cosé ) Pe sate, sing sind 
rcosd r 


1 
r 








o(- -2v,0, + Pe cos’ a sin & cosd fe » sind cos 


Thus, the Hamiltonian is 


H(x,u,A) =F +A'f « u) 
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H=1+2,7+A,0 +2,0 +4,, [6 cos’ + rb” — 4. fe cos ‘a 


Aus 


rcosd 





+ 


E 28 (1 sing — 00s 9) + Ecos? asinatsind | 


ae *(( ~2n + a cos’ a sin a cosd ] sing cose 
5 


The Lagrangian is therefore 
L(x(t), u(t), AQ) = (x(t), Wt), AD) +o, | (g(t), uD ) 
or simply, 
L=H +u,5 +U,0 where we 2.2] and [1s =0 if 5 is unconstrained. 
However in the actual implementation of control bounds, 6 is limited by 
[ x] 


de rz | to avoid non-distinct controls where two (0,5) control pairs could produce 


the same resulting thrust vector. It will be apparent later why this restriction is important. 


First deriving the control law for the clock angle we calculate 


ae he Bu cos’ a sin acos § — de But 
05 rcoso r? rr 





—cos’a sinasin§ +H, =0 


Pe 


Dividing by cosa sina , we obtain the tangent steering law for the clock 


angle 6 interior controls for «0 #0,4— 


uy 
t rey peepee aa = 
(29) an Xi cos > Hs 0 
The control angle 6 has no effect on the resulting thrust vector when the sail is 


exposing its full area to the sun or when it is exposing no area [ =0.4% Applying 
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the KKT conditions (Appendix A) to obtain the circumstances when the clock angle is at 


the “stops” we get 
T 
5 a ks. <0 
T 
5 = baal 20 


where the sign of the covector determines which way the clock angle is oriented. 


Next, we desire the steering law for the cone angle, 7. Applying the same KKT 


conditions to cone angle, , we obtain similar results. The “stops” of the cone angle 
control occur at = and . when the control dual variable meets the following 


conditions. 


Continuing with the Minimum Principle, we obtain 


ue =-i,, i = (3cos" 0. sin oc ) +—8— hg Be sind (-2cosa sin? O + cos *a) 
0a rcosd r 
os ” cos 6 (-2cosa sin? +cos a) +E, =0 
is 
Dividing by cos*a produces 





r 
(30) -3A,, tana + dro sind (—2tan? a +1)+—*cos8 (—2tan’a +1)+p, =0 
rcos > r 
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T ‘ Ste De 
where a #+ a This corresponds to the interior controls only when the cone angle 


multiplier is zero, 1, =0. Using equation (29) for tand, we form an expression for cos} 
and sin6 as a function of the states and costates. 


Ayo COSO 


cos6 =+ sin6 =+ 


yp COS "2 Avy COS eae, 
® ® 


The sign ambiguity is resolved when we take into account the control bounds, i.e. 


de -E.| , therefore sind > 0. Recall that equation (29) is not valid at a =0. 


Rearranging equation (30) into standard quadratic polynomial form for tana 
2 m 1 
atan°-a+btano+c=0, aE (-2.5] and a #0 
where the coefficients are given by the following: 


nr 
a=-2 ho sind —2—“cos8 
rcosd r 








) ae r 
c=—*® sind +—*cos 8 
rcos > r 


Solving the quadratic yields the tangent steering law for the cone angle interior 


controls. 
oa +4] 2 
(31) tang, = VST MC for ve(-2.5) and a #0 
a 


Following numerical solutions to simple optimal control problems (i.e. no interior 
knots), DIDO will output the time histories of costate as well as dual variables associated 


with the control angles. Steering laws in equations (29) and (31) use these outputs to 
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generate control profiles that may be compared to the DIDO control history output for 


verification. 


D. RESULTS 


1. Benchmark Problem Solutions with Elliptic Inclined Orbits 


The question in reviewing the same series of benchmark problems with elliptic 
inclined orbits is the following. Does adding a degree of freedom in the dynamics to 
meet the 3D event conditions enable us to further reduce the cost from the elliptic 
coplanar case or does it hurt us? Apparent from the first several solutions is that the 3 
DOF model minimizes the transit time to a smaller value than the circular coplanar 


model, but produces a larger trajectory time than the elliptic coplanar case. 


a. Earth-Mars Flyby 


To reach Mars as quickly as possible and sail past without gravitational 
interaction, it was intuitive that Mars periapsis provided the minimum path distance in the 
coplanar model. Now when we consider that the target orbital planes are inclined with 
respect to each other it may be desirable to choose a path that avoids orbit “cranking” 
even when hitting the periapsis is not possible in plane. The solution shows that the sail 
changes planes a small amount in order to reach Mars at its perihelion as shown in Figure 
44. Mars’ orbit has a small inclination with respect to the ecliptic so it is preferable to 
traverse a short distance even though a small cranking maneuver is incurred. Figure 45 
reveals that the time to intercept Mars with our standard sail ( B=.17) is 138 days, which is 
about the same as transiting to an in-plane elliptic Mars orbit (137 days). The sail cone 
angle control history is very similar to the coplanar case cone history while the clock 
angle steadily increases (Figure 46). Controls derived from the tangent steering law 


using costates are compared to the DIDO control output showing a close match. 
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Mars flyby sail trajectory. DIDO states and propagation 





Figure 44 —_—‘ Flyby Mission to Mars with Elliptic, Inclined Planetary Orbital 
Planes. DIDO output (dots) and propagated path (line). Mars orbit inclination is 
exaggerated for display purposes. 
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Mars flyby. DIDO states and propagated states. 20 nodes. 





L 
r= 
a 
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Zz 
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time acer 
Figure 45 State History for Mars Flyby (3D Model) 
Steering law with Costates vs. DIDO output. 
300 
250 ——-- tan steering cone 
——-- tan steering clock 
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DB 0 dido clock 
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Time (years) 


Figure 46 Cone and Clock Angle Controls for Mars Flyby. History is shown for 
DIDO (markers) and tangent steering law using states and costates (lines). 
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b. Rendezvous 
Observing the and @, states in Figure 47 and Figure 48 gives insight 


into what is required to perform the plane change maneuver. The sail trajectory remains 


in the Earth ecliptic plane for the first quarter of the mission time, then “cranks” the 


orbital plane to match the final sail and @, with Mars exactly at the time of Mars 


encounter. The mission is 1.01 years, which is a month faster that the circular coplanar 


and takes only marginally longer than the elliptic coplanar orbit model (0.977 years). 


80 nodes 


ND units 





“9 0.2 0.4 0.6 0.8 1 1.2 1.4 
Time (years) 


Figure 47 EM Rendezvous with Elliptic, Inclined Orbits. 
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80 nodes 


ND units 





Time [years] 


Figure 48 Optimal Profile for @ and @, States in the 3D Rendezvous Mission. 


Turning to the costate analysis, it is evident that as in the 2-D solutions 
there is a critical time during the maneuver at which all costates are at a local maximum 
or minimum or zero (Figure 49). It is beyond the scope of this thesis as to why this 
behavior occurs, however at this point in time (about 1/4" into the mission time) the 
control angles make large adjustments to achieve the state changes previously described. 
When inserted into the control laws derived in equations (29) and (31), the costates and 
states will produce derived cone and clock angle controls. A comparison of costate- 
derived controls and DIDO controls is shown in Figure 50. The optimal control does not 
rotate the cone angle to the extremes that the coplanar model control did. The controls 
here only pressed the cone angle between 18 and 60 degrees versus between 10 and 70 
degrees. The optimal cone angle has to account for the clock angle as well in this 3 DOF 


model. 
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Costates for the minimum time 3-D Mars rendezvous 






ND units 
= 
Oo 


-150 


-200 


Time [years] 


Costate History for Minimum Time Rendezvous with 3D Orbits 


Figure 49 
Model. 


Checking feasibility, the propagated controls were observed to match 


closely with the DIDO output (Figure 51). For optimality, the Hamiltonian is constant in 


Figure 52. 
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Derived controls using costates and DIDO controls 


rity 
Beer tiiiiey 


—_ 
ee 


control angle (deg) 





time [years] 


Figure 50 Control History for 3D Rendezvous. Displayed are DIDO controls 
(markers) and the controls derived from the tangent steering law with the states and 
costates as inputs (lines). 


Propagated path and DIDO state output 





Figure 51 Propagated Path of 3D Rendezvous (line) with DIDO Output (dots). 
Mars inclination is exaggerated for display purposes. 
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Hamiltonian for 3D EM rendezvous. 80 nodes 
aa 


2125: 


2 


-2.5 


3 


-3.5 


0 0.2 0.4 0.6 0.8 1 1.2 
Time (years) 


Figure 52. Hamiltonian of 3D EM Rendezvous Problem 


Cc. Double Rendezvous 


To optimize the Earth-Mars-Earth (EME) double rendezvous, we add to 
the rendezvous event conditions a return trip to Earth. The 2-D elliptic coplanar double - 
rendezvous exhibited symmetry in states and controls on the outbound and inbound legs 
(Figure 31). Now using inclined ellipses with their corresponding positions and velocity 
vector fields as target manifolds, we seek insights into the double -rendezvous with this 


higher fidelity model. 


The event conditions now include an interior event corresponding to a 
rendezvous with Mars and a final event corresponding to a rendezvous with Earth. To 
ensure Earth is actually present when the sail intersects its orbit manifold, we constrain 
the time of flight using Kepler’s equation. This last non-linear constraint replaces the 
corresponding one in the 2-D model that was listed in Table 4. Fortunately, the non- 
linear equation does not need to be solved explicitly as it simply serves as a constraint to 
the trajectory. In this way, the time of flight is coupled with the eccentric anom aly that 


the Earth traverses, E. Earth’s angular displacement traversed is constrained by the 
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flight time of the spacecraft. As long as Earth’s true anomaly change is the same 
displacement as the spacecraft angular displacement plus a multiple of 2% phasing will 


have been achieved. 


In order to rendezvous with a planet, we need to describe the position and 


velocity of both spacecraft and planet in 3-D (Np=3) for all 4 events (Np=4). 


Spacecraft 3D position at events in Eframe Cartesian coordinates are 


expressed as 


Ry =| Yo Wo Yi Ye 


and likewise spacecraft 3-D velocity at events in LVLH. 


- + 


V0 v ri Vii Vit 


= - + 
V.=! Veo Vi Yer Vo f 


Transforming the sail velocity vector from LVLH into in E-frame, we use 


the transformation matrix. 


B=[B] =[R,(-9)][R,@)] 


23 
So spacecraft velocity in E-frame Cartesian coordinates is 
= NpXNg 
Ve = BY, Vie eR ‘ 


Because the craft is returning to Earth, it must be properly phased with 
Earth to ensure the planet is there when the spacecraft arrives. This constraint makes an 
Earth orbital motion model necessary accounting for its position in space and time. The 
following equations are used to constrain initial and final Earth events to ellipses in three 


dimensions. All state vector multiplication and division operators are element -wise. 


First, we define the Earth initial and final angular positions in the perifocal 


frame. 
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0.=[0.. O, |=[0, 6, |-o.[n, 2], 


Notice that 0, ¢ R* where [9, 0 | are the initial and final angular 


spacecraft positions in E-frame which are coincident with Earth’s initial and final 
position. The following equations relate Earth’s eccentric anomaly to its true anomaly. 
Eccentric anomaly is expressed as a vec tor with each element corresponding to a knot, or 
end condition. The purpose of expressing it in this manner is to be able to handle 


multiple cycles with multiple encounters ( NV, >1) using a single variable. 


The initial and final cosine of Earth’s eccentric anomaly are contained in 
the vector cosE. The vector operations are performed in an “element-wise” fashion 


where the cosine of the vector ®, is the cosine of each individual element of vector 8. 


e, +cos (0, ) ‘ 
cosE = ————~_,, cosEe R™ 
1+e, cos (8, ) 


Earth’s radial position as a function of cosE is given by 
r,=a,(l-e,cosE) ,r,¢ R™ 
where r, =[r,).7,,]' so that vector sinE may be expressed as 


_ r,sin(O,) 
a,Jl-e,’ 


Thus, Earth’s initial and final eccentric anomaly are contained in the 


sinE , re R™ 


vector E. 





E= tan" (= ) , Ee R” Earth’s eccentric anomaly 
cos 


The next set of equations is used to define Earth’s position based on the 


initial and final time. 


n= |= Earth’s mean motion 
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M,=£,+e,sinE, Initial mean anomaly 


M , =" bopansir Final mean anomaly 


Now we define the initial and final Earth position event matrix in 


heliocentric Cartesian coordinates, E-frame where. 


r, cos(8,) 


cos(0 p 


i 
Rx =| 7, sin®,) 7,sin@,) 





The speed of Earth on its elliptical path around the sun is 











where cos, and sinB, are the cosine and sine of Earth’s flight path angle at the event 
knots. Earth’s velocity in radial and transverse components in the perifocal frame are 
expressed in the following equations where the symbol ® is used to emphasize that this 
is an “element-wise” multiplication such that each element in one vector is multiplied by 


the corresponding element in the other vector. 














u, =V, @sinB, ,u,e R” 
v, =V, @cosB, ,v,e R™ 
u, 
Neer = [ R,(-8, ) Ms > Vucarr € Ries 
0 
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= NpXNez 
Vex = AY cane » Vee RO 


where A, is a matrix performing a 3-1-3 transformation from Earth perifocal to Cartesian 
heliocentric inertial (E-frame). See Appendix D for the derivation. 
Having defined all the necessary quantities, it is now a simple matter to set 


boundary conditions. Letting the spacecraft initial and final event states be Rgzna and 


VsEnd We have 


f X) x f Y. Xp 
Rona Yo ye > Viena Vs Vp 
% V, V:, 


[° a i a 
R,=l\y yl. Vay Vy 


G Zr Vv OV 


Ry —R « =9 
and 
Views a =0 
We can establish the Mars intermediate events in a similar fashion. 
R,-R,,, =0 
and 
V,-V,,=0 


One last event constraint ensures proper phasing with Earth. This 3s 


accomplished by forcing Kepler’s equation as a constraint. 
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E-e,sinE-M, = =0 


Earth mean anomaly and eccentric anomaly are represented by M and E 
respectively. The event constraints for the 3D EME rendezvous mission are summarized 


in the table below. 


Double rendezvous 
constraints 
[lower,upper] 


[RyRy] 


ce? 


[V 


eE? 


Vir ] 


Table 11 EME 3D Double Rendezvous Event Constraints 





States and controls for the 3 DOF model are bounded by 


xeScR* 


ucUcR 


Where the set S is chosen to avoid singularities in the dynamics and the set 
U was chosen to avoid duplicity in the controls. Because it was desired to obtain distinct 


(a 6 ) pairs (every sail normal vector orientation is represented by only one (6 ) pair), 
the bounds were such that 0 <a <t and -m <5 <z. It was recognized that the sin and 


cos components of the normal vector could have been used as the controls with proper 
path constraints, however the numerical solutions took longer to complete, so the cone 


and clock angles were retained as the controls. 
Results and Analysis 


Once more, Mars aphelion played the role of the optimal rendezvous point 


at the intermediate event for the round trip. Modeling Mars’ small inclination with the 
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ecliptic produced a solution that looked very similar to the 2D ellip tic when case when 
looking from the top view (compare Figure 53 with Figure 31). The mission time, 
however, was slightly longer than the 2D case taking 2.42 years instead of 2.3 years since 
the spacecraft must change planes twice. The mission time of the 3D model is the same 
as that of the 2D circular orbits model except that there is no longer the precise symmetry 
of states and controls exhibited by the other model. In three dimensions, the spacecraft 
takes 1.03 years on the outbound leg and 1.39 on the inbound (2D circular coplanar 
model required 1.205 years/leg). One of the amazing features of the path is that the 
spacecraft sails out of both Mars and Earth orbital planes for a significant portion of t he 
return trip to Earth rendezvous (igure 54). This large departure from the planetary 
planes is accomplished in order to reach Earth orbit in phase with Earth to an accuracy 
2.64E-4 radians (~39E3 km from the center of the Earth). Testing the feasibility, the 


control profile was entered into the propagator producing the result in Figure 55. 





Figure 53 3D EME Double Rendezvous Path (top view) 
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Figure 54 3D EME Double Rendezvous Path (oblique view) 


Figure 55 Propagated Path (line) and DIDO States (dots) for 3D EME Double 
Rendezvous 
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2. Earth-Mars Cycler Set-up 


Formulating the 3-D elliptic orbit cycler model into an OCP is a formidable task 
because cyclic end conditions are no longer imposed every EME round trip. Many round 
trips occur before Mars and Earth positions are exactly the same as they were at the start 
of the trajectory. The events as structured in the previous section are set up to handle 
multiple passes for a “cycle”. Since meeting synodic conditions with each round trip is 
no longer required, a possible OCP formulation would be to seek the minimum time to 
achieve a given number of EM passes. Using the framework outlined in this thesis, this 
would entail performing successive solutions to increasing number of EM visits (and thus 
number of knots) using intercept end conditions until initial state conditions are achieved 


at the final time. 
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Vv. CONCLUSIONS AND FUTURE WORK 


A. CONCLUSIONS 


Solar sails provide a viable means of propulsion for an Earth-Mars cycler. For a 


solar sail with lightness number B =.17, an orbit can be maintained with rapid revisit 
times between destination planets and slow Vs at each planet, both desirable cycler 


characteristics. Furthermore, the cycler lifetime is only dependent on sail degradation in 
the space environment since no propellant is depleted making solar sail cyclers attractive 


in comparison to conventional low thrust and chemical propellant cyclers. 


The framework in this thesis outlines a method of state and control optimization 
for both 2-D and 3-D cycler models. This framework can be applied to other missions as 
well by changing the event conditions. With relatively simple dynamics and sail models, 
trajectory design could be accomplished by properly formulating desired events in the 
form of constraints. The pseudospectral method used within DIDO to solve the OCPs 
produces optimal controls and states that were verified as feasible (through propagation 
of initial conditions with a third party ODE solver) and optimal using the necessary 


conditions for optimality. 


B. OPPORTUNITIES FOR FUTURE WORK 


The following are some areas of interest for continuing work. 


Earth-Venus-Mars cycler. Is a cycler possible with three target planets and if 
so, what does it look like? Is it possible to reduce the revisit times from a two-planet 
cycler? This problem would require the solution to include the optimal order of planetary 


encounters. 


Use a solar sail to establish an optimal halo orbit. Many non-Keplarian orbits 
are possible with solar sails. Orbits around Lagrange points may be modeled and 
examined for missions such as early detection of geomagnetic storms on the sunward side 
of the Earth-Sun L1 point (ref 9 p. 231). Cycling orbits between L points could also be 


explored. 
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Implement two synodic period cyclers with solar sails. Reference 19 outlines a 
plan to make multiple passes of a target planet within a cycle (2 synodic periods). Planet - 


centered dynamics may play a role if trajectories spend significant time in a planet’s SOI. 


SEP/Solar Sail hybrid cyclers. Modeling a SEP/solar sail hybrid spacecraft 
would show when it is optimal to favor the SEP and when it is preferable to favor the sail 


for various cycler orbits (or any trajectory for that matter). 


Variable Isp cycler. When is it optimal to use high thrust-low Isp and when is it 
best to use low thrust-high Isp in a cycler trajectory? There are trades to examine with 


weighted cost functions between fuel mass and time. 


Make smoother controls near hard knots. Changing the controls to include 
“inertia” assisted in making controls somewhat smoother near knots, however for 3D 
problems it was important to implement constraints to ensure that arbitrary sail 
orientations resulted from unique cone angle/clock angle pairs. Implementing the sin and 
cos components of these angles as controls (with a path constraint) slowed the algorithm 
down and was not implemented, however this may have alleviated some of the jumpy 


control outputs. 


Improve gravity assist model. The matched asymptotes model as presented in 
this thesis is accurate for the problems posed, however more accurate solutions to some 
problems may call for a better gravity assist model. Possibilities are to either change 
dynamics inside SOI or use a jump in time and theta at the in terior knot in addition to 


velocity. 


Improve sail model. A higher fidelity model of the sail is available that includes 
billowing and diffuse reflection (see ref 9 pp.51-54). Once the dynamics and events 
models have been thoroughly tested, there are some design trades that could be explored 


as far as the sail itself. 


Change to dynamic equations that are better suited for cyclers. The 
discussion on coordinate systems outlined some benefits and pitfalls of using certain 
coordinates in cycler trajectories. Perhaps other coordinates than those presented here 


can be better implemented for optimal cyclers (e.g. equinoctial). 
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Solar sail stability. Solar sails have been used in spacecraft to dump 
accumulated momentum in reaction wheels due to torque perturbations. A possible OCP 
would be to minimize the time to dump momentum in reaction wheels for a given earth - 


centered orbit. 
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APPENDIX A. APPLICATION OF THE MINIMUM PRINCIPLE??” 


The first step in solving an optimal control problem is to construct the scalar 


Hamiltonian function, H , 
H(x,u,t,A) =F (x,u,t)+A(t) f (x,u,7) 


where f(x,u,t) are the dynamic constraints on the system, and A(t) are the Lagrange 


multipliers called costates. According to Pontryagin’s Minimum Principle, the optimal 


control is obtained by solving the following problem at each instant in time. 


Minimize H with respect to u, with ueU where U is the set of all allowable 


controls. To solve this problem, the Lagrangian of the Hamiltonian is written as: 
L(x,ut.A.M, ) =H (x,u t,A)+M, (t)' g(x,u,t) 


where g(x,u,t) are the path constraints and p, (r) are the associated path covectors. 


The path constraints include all trajectory path constraints as well as any state and control 
bounds. Applying the Karush-Kuhn-Tucker (KKT) theorem to the Lagrangian results in 
a set of first order necessary conditions and provides a means to demonstrate the 


optimality of a solution. 


(32) oe =0 


with 
<0 8, = g(x,u,T) 
20 g(x,u,T)= g, 
H. = if ( ) 
=0 8, <8 (X,U,T)<g, 
any §1= 8, 
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The third case above describes a special condition when the constraints in g are 


interior or non-binding, 
8, <g(x.u.t)<g, 
For these cases, the multipliers, py, =0 and equation (32) simplifies to 
OE 20H 2 
du du 


A necessary second order condition for optimality when p, =0 is that 


(33) 
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APPENDIX B. DERIVATION OF THE 2-D EQUATIONS OF 
MOTION 


Lagrangian mechanics are used to achieve the 2D equations of motion using 
polar coordinates. Velocity variables are then written into a more familiar form using the 
radial and transverse velocity components. Figure 1 shows the states and coordinate 


convention used in the derivation. In general form Lagrange’s equation is written as 


where L 8 the Lagrangian, x is the generalized coordinate vector, and Q). is the 
generalized non-conservative force vector'®. 
We start with writing the specific kinetic energy as 


T =ty-yah or? +16?) 
2 2 


Potential energy due to gravity is 


The Lagrangian may be written as 


eT evel ga reryst 
2 r 


The non-conservative generalized force is due to the solar radiation pressure 


(SRP) on the solar sail. Acceleration magnitude is given by 


T_ Bu 


—_—= cos’ a 
m % 


Thus the acceleration in the radial direction is 


T 
a, =—cosa = BU att 


r 


m % 


and the applied acceleration in the along-track direction is 
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Le : 
dy =—sina = BU eG sina 
m i) 


Writing the Lagrangian equations of motion, we obtain 


d( ol = OL _ ig? _ th 
dt\ or or r 


(34) ¢—07+ # 


im= 
r r 


= oe =2ri0+r0 ae, 
dt\ 08 ed 


(35) 2rr0 +78 =a,r 


where a,r is a generalized specific torque. 


Equation (35) reduces to the conservation of angular momentum when a, =0. 


Now to convert the 7s and @s into u and v velocity components, we use the 


following relations. 

Transverse velocity is 

(36) v=r0 
and transverse acceleration is 

(37) b= 10 +70 


Using the relation in equation (36) with equation (34) we write the equation of 


motion for . 


(38) 





Rewriting equation (35) and factoring an r yields 
(39) r(r6 +70 +70 )=agr 
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Recall that u=7 therefore w = 770 . Using (37) with (39) we obtain the equation 


of motion for v. 


rv +r1® = a,r 


uw 
b=-—+a, 
ow 
(40) v= SM oe rey 
rt 
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APPENDIX C. DERIVATION OF THE 3-D EQUATIONS OF 
MOTION 


As with the 2-D case, we use the Lagrangian approach to acquire the 3-D equations of 


motion. 


Figure 56 shows the RSW spherical coordinate system (ref 4 p. 42) and 


convention used to represent the states. 





Figure 56 © RSW Coordinate System 
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In general form Lagrange’s equation is written as 
(41) ets r 


where L is the Lagrangian, x is the generalized coordinate vector, and Q,. is the 
generalized non-conservative force vector'®. The generalized non-conservative force for 
this case is the solar radiation pressure (SRP) on the solar sail resolved in radial (R ), 


along-track (S ), and cross-track (W ) directions. 


We start by writing the specific kinetic energy as 
r = iy “Vv 
2 


where the velocity v is 
v=iR+0 cos oS + roW 


Specific potential energy is 


making the Lagrangian 
17. : 2 +\2 wu 
L=T-V=-|i+ r8cosd) +(7 | 
1724+ (rcoso)' +(r)'] +! 


The external acceleration due to the SRP on the sail has a magnitude of 


_T _ By 


a=—=-~cos’o 
m or 


and have components in the radial, along-track and cross-track directions, written as 


acosaR+asind sin S+asina cosé W 


Having established the necessary quantities, we proceed to derive the Lagrangian 


equations of motion. For the r state, we have the following. 


d oL a OL ws) 9 <2, u 
—_|] — ps —=/70 _— 
naka fi or eat r° 
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(42) # — 10° cos’ d—ro7?+ P= acoso 


r 
Turning to the @ state we form the quantities 


alas 70 (-2 cos sind )+cos” (r6 +2770 ) cae, 


The generalized force is a torque about the W axis 


Q,, =asina sin 8 (r cos >) 
Lagrange’s equation (41) becomes 
re (-2p cos sing }+cos’o (r°6 +2ri0 ) = asino sind (rcos¢ ) 
Dividing by rcoso@ results in 
) (—2sinod )+ coso7@ + 278 = asina sind 


Finally, we rearrange terms to get the equation of motion for 6 . 


(43) 


[ 20 (1 sind — rcosd ) +asinasind | 


rcoso 





Note that we need to be cautious at the singularities where r=0 and cos o=0. 


Now for the @ state we write 


d(oL se bes OL as . 
ose err 3p r0* cosd sing 


The generalized force in the @ sense direction (torque in this case) is given by 
Q, = arsina cosd 
Substituting into equation (41) yields 


2rrd +r°o +r? cos sind = ar sinacos 
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Dividing by r and rearranging gives the equation of motion for o ; 


(44) o = Ly 27 + asina cosd ) —0’ sind cosd 
r 
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APPENDIX D. TRANSFORMATION INTO HELIOCENTRIC - 
ECLIPTIC COORDINATES (E-FRAME) 


The matrix derivations below are used to transform the spacecraft state and event 
conditions into a common coordinate system — heliocentric ecliptic Cartesian coordinates, 
called the E-frame. Spacecraft states at the events (position and velocity) are given in 
spherical coordinates while the target planet manifolds (elliptical paths with velocity 


related to position) are most easily expressed in a perifocal frame. 





x,y 


Figure 57 Spherical coordinate system 


SPACECRAFT 


Referring to the Figure 57, we transform the spherical position coordinates into 
Cartesian coordinates as follows. It is worth noting that some texts (e.g. ref 18) choose 
the @ coordinate as the complement angle from the one depicted above, so it is important 


to be careful in deriving the transformation matrices. 


rcosd cos® 


rcosd sin® 
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The spacecraft velocity state has elements given as 7,0,0. In the parlance of 
DIDO, these spacecraft states are the “primal.states”. First we need the velocity vector of 


the spacecraft in v,,%4,v, components where the subscripts indicate which sense 


() 


direction the velocity vectors point. This transformation appears as 


lv, | i 
v, | =| 18 cos 
Vy ro 
or in our notation 
V, V, 
(45) Vy |=| 7, Cos 
V% 1, 


These coordinates may now be transformed into the E-frame using a 3-rotation 


and 2-rotation matrix. 


»,) =[R,(9)][R, 6)]] 


We choose to call the matrix | R, (6) |[ R,(o) | the B matrix. 


cos®8 cos -—sin@ cosOsingd 
sind cosd = cos8-—s sin sind 


—sing 0 cos 





Using transformations given in boxes, we can transform the position and velocity 
of the spacecraft in spherical coordinates (primal.states) into the Eframe coordinate 


frame. 
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PLANETS 


To define the manifolds of Earth and Mars we need their positions, defining an 
ellipse, and equations for their velocities at any point on the ellipse. Time is imposed as a 
non-linear constraint in the event conditions without need for transformation, so it will 
not be addressed here. An ellipse may be represented easily in either polar or Cartesian 
form, so we choose the polar form since most astrodynamic textbooks use this form. 
Positions and velocities in the perifocal frame (Figure 58) are related by the following 
equations. 


r=—P and V= wu oe 
1+ecos0 , 


The 6, used in these equations is the angular positions of the planet from their 


perihelion at a particular event. The relation between the planet true anomaly and the 


spacecraft angular position is 


9, =8 -(@, +Q,) 


Pp 
for small inclinations (ref 4). For larger inclinations another transformation is required. 


Taking the locus of planet positions, the ellipse, and transforming in to Cartesian 


coordinates we write 


Xx, =rcos6, 
y,=rsine , 


z, =0 


Pp 


Transforming into the E-frame requires a 3-1-3 rotation matrix defined using the 


longitude of ascending node Q, inclination i, and argument of periapsis @ (ref 13). 


This matrix that transforms Cartesian coordinates from the perifocal frame to the 


E-frame is defined as 


A =[R,(-o)][R Ci) ][R (-2)] 


Eep 
cosQcos@—cosisinQsin@ —-cosisinQcos@—cosQsin@ — sinisin@ 
=| cosicosQsin@+sinQcos@ cosicosQcos@—sinQsin@ —sin icosQ 
sin isin @ sin icos@ cos i 
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R denotes rotation matrices about the 1, 2, or 3 axis. See ref 13, p.80, ref 18, 


p.190, and ref 4, p. 151 for more details. 


Now transforming the planetary velocity vectors from the perifoc al coordinates 
(Figure 58), u and v in figure to the E-frame. First, transform u and v into Cartesian 


coordinates in the perifocal frame. 


V,. u 
Vv, =[R,(-9,) v 
Vv, perifocal 


Finally, transforming the Cartesian velocity vector in the per ifocal plane into the 


E-frame, we use the A matrix again. 


Vy Vy 
vy} =Ayy, 
Vv, E V, perifocal 





Figure 58 Planetary perifocal system 
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APPENDIX E. SOLAR SAIL CONTROL HODOGRAPH ANALYSIS 


A hodograph is useful in determining convergence char acteristics of a solution 
using a particular set of controls in the dynamics. Figure 59 shows the hodograph, the 
mapping of the control, @, onto the w—v plane. Because the thrust magnitude and 
resulting acceleration due to the sail is completely coupled with the sail angle, the 
controls map as closed curve, not a region. Equations (38) and (40) plotted for all 
feasible & values produces the hodograph. The curve is not convex in that a line drawn 
between any two points on the curve does not remain completely in the locus of points 
defined by the control mapping. There may be a relation of convexity and convergence 


speed. 


——= 


Figure 59 Hodograph for Cone Angle, 


Using a new control scheme by setting the control variable equal to & maps as a 
straight line, which is inherently convex using the above definition. Solutions do in fact 


converge faster for this case. 
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